Belle II Software development
SVDdEdxValidationAlgorithm.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9#include <svd/calibration/SVDdEdxValidationAlgorithm.h>
10
11#include <tuple>
12#include <vector>
13#include <string>
14
15#include <TROOT.h>
16#include <TStyle.h>
17#include <TMath.h>
18#include <TFile.h>
19#include <TColor.h>
20#include <TLegend.h>
21#include <TCanvas.h>
22#include <TH1D.h>
23#include <TH1F.h>
24#include <TH2F.h>
25#include <TAxis.h>
26#include <TGraph.h>
27#include <TMultiGraph.h>
28#include <TCut.h>
29
30#include <RooDataSet.h>
31#include <RooRealVar.h>
32#include <RooAddPdf.h>
33#include <RooGaussian.h>
34#include <RooChebychev.h>
35#include <RooBifurGauss.h>
36#include <RooDstD0BG.h>
37#include <RooAbsDataStore.h>
38#include <RooTreeDataStore.h>
39#include <RooMsgService.h>
40#include <RooStats/SPlot.h>
41#include <ROOT/RDataFrame.hxx>
42
43using namespace ROOT;
44using namespace RooFit;
45using namespace Belle2;
46
47SVDdEdxValidationAlgorithm::SVDdEdxValidationAlgorithm() : CalibrationAlgorithm("SVDdEdxValidationCollector"),
48 m_isMakePlots(true)
49{
50 setDescription("SVD dE/dx validation algorithm");
51}
52
53/* Main calibration method */
54CalibrationAlgorithm::EResult SVDdEdxValidationAlgorithm::calibrate()
55{
56 gROOT->SetBatch(true);
57
58 // const auto exprun = getRunList()[0];
59 // B2INFO("ExpRun used for calibration: " << exprun.first << " " << exprun.second);
60
61 // Get data objects
62 auto TTreeLambda = getObjectPtr<TTree>("Lambda");
63 auto TTreeDstar = getObjectPtr<TTree>("Dstar");
64 auto TTreeGamma = getObjectPtr<TTree>("Gamma");
65
66 if (TTreeLambda->GetEntries() < m_MinEvtsPerTree) {
67 B2WARNING("Not enough data for calibration.");
68 return c_NotEnoughData;
69 }
70
71 // call the calibration functions
72 TTree* TTreeLambdaSW = LambdaMassFit(TTreeLambda);
73 TTree* TTreeDstarSW = DstarMassFit(TTreeDstar);
74 TTree* TTreeGammaWrap = TTreeGamma.get();
75
76 std::vector<TString> PIDDetectors;
77 PIDDetectors.push_back("SVDonly");
78 if (m_fullValidation) {
79 PIDDetectors.push_back("ALL");
80 PIDDetectors.push_back("noSVD");
81 }
82
83 std::map<TTree*, TString> SWeightNameMap = {
84 {TTreeGammaWrap, "1"},
85 {TTreeDstarSW, "nSignalDstar_sw"},
86 {TTreeLambdaSW, "nSignalLambda_sw"}
87 };
88
89 for (const TString& PIDDetectorsName : PIDDetectors) {
90 PlotEfficiencyPlots(PIDDetectorsName, TTreeGammaWrap, SWeightNameMap[TTreeGammaWrap], "FirstElectron", "electron", TTreeDstarSW,
91 SWeightNameMap[TTreeDstarSW], "PionD", "pion",
92 "BinaryElectronPionID",
93 "0.5", m_NumEffBins, 0., m_MomHighEff);
94 PlotEfficiencyPlots(PIDDetectorsName, TTreeGammaWrap, SWeightNameMap[TTreeGammaWrap], "FirstElectron", "electron", TTreeDstarSW,
95 SWeightNameMap[TTreeDstarSW], "Kaon", "kaon",
96 "BinaryElectronKaonID", "0.5",
97 m_NumEffBins, 0., m_MomHighEff);
98 PlotEfficiencyPlots(PIDDetectorsName, TTreeLambdaSW, SWeightNameMap[TTreeLambdaSW], "Proton", "proton", TTreeDstarSW,
99 SWeightNameMap[TTreeDstarSW], "PionD", "pion",
100 "BinaryProtonPionID", "0.5",
101 m_NumEffBins, 0.25, m_MomHighEff);
102 PlotEfficiencyPlots(PIDDetectorsName, TTreeLambdaSW, SWeightNameMap[TTreeLambdaSW], "Proton", "proton", TTreeDstarSW,
103 SWeightNameMap[TTreeDstarSW], "Kaon", "kaon",
104 "BinaryProtonKaonID", "0.5",
105 m_NumEffBins, 0.25, m_MomHighEff);
106 PlotEfficiencyPlots(PIDDetectorsName, TTreeDstarSW, SWeightNameMap[TTreeDstarSW], "PionD", "pion", TTreeDstarSW,
107 SWeightNameMap[TTreeDstarSW],
108 "Kaon", "kaon",
109 "BinaryPionKaonID", "0.5", m_NumEffBins,
110 0., m_MomHighEff);
111 PlotEfficiencyPlots(PIDDetectorsName, TTreeDstarSW, SWeightNameMap[TTreeDstarSW], "Kaon", "kaon", TTreeDstarSW,
112 SWeightNameMap[TTreeDstarSW],
113 "PionD", "pion",
114 "BinaryKaonPionID", "0.5", m_NumEffBins,
115 0., m_MomHighEff);
116 }
117
118 if (m_fullValidation) {
119 PlotROCCurve(TTreeGammaWrap, SWeightNameMap[TTreeGammaWrap], "FirstElectron", "electron", TTreeDstarSW,
120 SWeightNameMap[TTreeDstarSW], "PionD",
121 "pion", "BinaryElectronPionID");
122 PlotROCCurve(TTreeGammaWrap, SWeightNameMap[TTreeGammaWrap], "FirstElectron", "electron", TTreeDstarSW,
123 SWeightNameMap[TTreeDstarSW], "Kaon",
124 "kaon", "BinaryElectronKaonID");
125 PlotROCCurve(TTreeLambdaSW, SWeightNameMap[TTreeLambdaSW], "Proton", "proton", TTreeDstarSW, SWeightNameMap[TTreeDstarSW], "PionD",
126 "pion",
127 "BinaryProtonPionID");
128 PlotROCCurve(TTreeLambdaSW, SWeightNameMap[TTreeLambdaSW], "Proton", "proton", TTreeDstarSW, SWeightNameMap[TTreeDstarSW], "Kaon",
129 "kaon",
130 "BinaryProtonKaonID");
131 PlotROCCurve(TTreeDstarSW, SWeightNameMap[TTreeDstarSW], "PionD", "pion", TTreeDstarSW, SWeightNameMap[TTreeDstarSW], "Kaon",
132 "kaon",
133 "BinaryPionKaonID");
134 PlotROCCurve(TTreeDstarSW, SWeightNameMap[TTreeDstarSW], "Kaon", "kaon", TTreeDstarSW, SWeightNameMap[TTreeDstarSW], "PionD",
135 "pion",
136 "BinaryKaonPionID");
137 }
138 B2INFO("SVD dE/dx validation done!");
139
140 return c_OK;
141}
142
143// generic efficiency and fake rate
144void SVDdEdxValidationAlgorithm::PlotEfficiencyPlots(const TString& PIDDetectorsName, TTree* SignalTree, TString SignalWeightName,
145 TString SignalVarName, TString SignalVarNameFull, TTree* FakeTree, TString FakeWeightName, TString FakeVarName,
146 TString FakeVarNameFull, TString PIDVarName, TString PIDCut, unsigned int nbins, double MomLow, double MomHigh)
147{
148
149 if ((SignalTree == nullptr) || (FakeTree == nullptr)) {
150 B2FATAL("Invalid dataset, stopping here");
151 }
152
153 if ((SignalTree->GetEntries() == 0) || (FakeTree->GetEntries() == 0)) {
154 B2FATAL("The dataset is empty, stopping here");
155 }
156
157 if ((SignalTree->GetBranch(Form("%sMomentum", SignalVarName.Data())) == nullptr)
158 || (FakeTree->GetBranch(Form("%sMomentum", FakeVarName.Data())) == nullptr)) {
159 B2FATAL("Check the provided branch name, stopping here");
160 }
161
162 TString SignalFiducialCut = "(1>0)"; // placeholder for a possible sanity cut
163 TString FakesFiducialCut = "(1>0)";
164
165 // Produce the plots of the SVD PID distribution
166 if (PIDDetectorsName == "SVDonly") {
167 SignalTree->Draw(Form("%s%s%s>>hSignalPIDDistribution(100,0.,1.)", SignalVarName.Data(), PIDVarName.Data(),
168 PIDDetectorsName.Data()),
169 SignalWeightName + Form("* (%sMomentum>%f && %sMomentum<%f)", SignalVarName.Data(), MomLow, SignalVarName.Data(), MomHigh), "goff");
170 TH1F* hSignalPIDDistribution = (TH1F*)gDirectory->Get("hSignalPIDDistribution");
171 hSignalPIDDistribution->Scale(1. / hSignalPIDDistribution->Integral());
172 hSignalPIDDistribution->GetXaxis()->SetTitle(PIDVarName + PIDDetectorsName + " for " + SignalVarNameFull);
173 hSignalPIDDistribution->GetYaxis()->SetTitle("Candidates, normalised");
174 hSignalPIDDistribution->SetMaximum(1.35 * hSignalPIDDistribution->GetMaximum());
175
176 SignalTree->Draw(Form("%sElectronLLSVDonly>>hSignalElectronLLDistribution(100,-17.,3.)", SignalVarName.Data()),
177 SignalWeightName + Form("* (%sMomentum>%f && %sMomentum<%f)", SignalVarName.Data(), MomLow, SignalVarName.Data(), MomHigh), "goff");
178 TH1F* hSignalElectronLLDistribution = (TH1F*)gDirectory->Get("hSignalElectronLLDistribution");
179 SignalTree->Draw(Form("%sPionLLSVDonly>>hSignalPionLLDistribution(100,-17.,3.)", SignalVarName.Data()),
180 SignalWeightName + Form("* (%sMomentum>%f && %sMomentum<%f)", SignalVarName.Data(), MomLow, SignalVarName.Data(), MomHigh), "goff");
181 TH1F* hSignalPionLLDistribution = (TH1F*)gDirectory->Get("hSignalPionLLDistribution");
182 SignalTree->Draw(Form("%sKaonLLSVDonly>>hSignalKaonLLDistribution(100,-17.,3.)", SignalVarName.Data()),
183 SignalWeightName + Form("* (%sMomentum>%f && %sMomentum<%f)", SignalVarName.Data(), MomLow, SignalVarName.Data(), MomHigh), "goff");
184 TH1F* hSignalKaonLLDistribution = (TH1F*)gDirectory->Get("hSignalKaonLLDistribution");
185 SignalTree->Draw(Form("%sProtonLLSVDonly>>hSignalProtonLLDistribution(100,-17.,3.)", SignalVarName.Data()),
186 SignalWeightName + Form("* (%sMomentum>%f && %sMomentum<%f)", SignalVarName.Data(), MomLow, SignalVarName.Data(), MomHigh), "goff");
187 TH1F* hSignalProtonLLDistribution = (TH1F*)gDirectory->Get("hSignalProtonLLDistribution");
188
189 // same but only for tracks that are expected to actually have SVD dEdx info
190 SignalTree->Draw(Form("%sElectronLLSVDonly>>hSignalElectronLLDistributionGood(100,-17.,3.)", SignalVarName.Data()),
191 SignalWeightName + Form("* (%sSVDdEdx>0) * (%sMomentum>%f && %sMomentum<%f)", SignalVarName.Data(), SignalVarName.Data(), MomLow,
192 SignalVarName.Data(), MomHigh), "goff");
193 TH1F* hSignalElectronLLDistributionGood = (TH1F*)gDirectory->Get("hSignalElectronLLDistributionGood");
194 SignalTree->Draw(Form("%sPionLLSVDonly>>hSignalPionLLDistributionGood(100,-17.,3.)", SignalVarName.Data()),
195 SignalWeightName + Form("* (%sSVDdEdx>0) * (%sMomentum>%f && %sMomentum<%f)", SignalVarName.Data(), SignalVarName.Data(), MomLow,
196 SignalVarName.Data(), MomHigh), "goff");
197 TH1F* hSignalPionLLDistributionGood = (TH1F*)gDirectory->Get("hSignalPionLLDistributionGood");
198 SignalTree->Draw(Form("%sKaonLLSVDonly>>hSignalKaonLLDistributionGood(100,-17.,3.)", SignalVarName.Data()),
199 SignalWeightName + Form("* (%sSVDdEdx>0) * (%sMomentum>%f && %sMomentum<%f)", SignalVarName.Data(), SignalVarName.Data(), MomLow,
200 SignalVarName.Data(), MomHigh), "goff");
201 TH1F* hSignalKaonLLDistributionGood = (TH1F*)gDirectory->Get("hSignalKaonLLDistributionGood");
202 SignalTree->Draw(Form("%sProtonLLSVDonly>>hSignalProtonLLDistributionGood(100,-17.,3.)", SignalVarName.Data()),
203 SignalWeightName + Form("* (%sSVDdEdx>0) * (%sMomentum>%f && %sMomentum<%f)", SignalVarName.Data(), SignalVarName.Data(), MomLow,
204 SignalVarName.Data(), MomHigh), "goff");
205 TH1F* hSignalProtonLLDistributionGood = (TH1F*)gDirectory->Get("hSignalProtonLLDistributionGood");
206
207
208 hSignalElectronLLDistribution->Scale(1. / hSignalElectronLLDistribution->Integral());
209 hSignalPionLLDistribution->Scale(1. / hSignalPionLLDistribution->Integral());
210 hSignalKaonLLDistribution->Scale(1. / hSignalKaonLLDistribution->Integral());
211 hSignalProtonLLDistribution->Scale(1. / hSignalProtonLLDistribution->Integral());
212
213 hSignalElectronLLDistributionGood->Scale(1. / hSignalElectronLLDistributionGood->Integral());
214 hSignalPionLLDistributionGood->Scale(1. / hSignalPionLLDistributionGood->Integral());
215 hSignalKaonLLDistributionGood->Scale(1. / hSignalKaonLLDistributionGood->Integral());
216 hSignalProtonLLDistributionGood->Scale(1. / hSignalProtonLLDistributionGood->Integral());
217
218 hSignalElectronLLDistribution->GetXaxis()->SetTitle(PIDVarName + "ElectronLLSVDonly");
219 hSignalElectronLLDistribution->GetYaxis()->SetTitle("Candidates, normalised");
220 hSignalElectronLLDistribution->SetMaximum(1.35 * hSignalElectronLLDistribution->GetMaximum());
221
222 hSignalPionLLDistribution->GetXaxis()->SetTitle(PIDVarName + "PionLLSVDonly");
223 hSignalPionLLDistribution->GetYaxis()->SetTitle("Candidates, normalised");
224 hSignalPionLLDistribution->SetMaximum(1.35 * hSignalPionLLDistribution->GetMaximum());
225
226 hSignalKaonLLDistribution->GetXaxis()->SetTitle(PIDVarName + "KaonLLSVDonly");
227 hSignalKaonLLDistribution->GetYaxis()->SetTitle("Candidates, normalised");
228 hSignalKaonLLDistribution->SetMaximum(1.35 * hSignalKaonLLDistribution->GetMaximum());
229
230 hSignalProtonLLDistribution->GetXaxis()->SetTitle(PIDVarName + "ProtonLLSVDonly");
231 hSignalProtonLLDistribution->GetYaxis()->SetTitle("Candidates, normalised");
232 hSignalProtonLLDistribution->SetMaximum(1.35 * hSignalProtonLLDistribution->GetMaximum());
233
234 hSignalElectronLLDistributionGood->GetXaxis()->SetTitle(PIDVarName + "ElectronLLSVDonly");
235 hSignalElectronLLDistributionGood->GetYaxis()->SetTitle("Candidates, normalised");
236 hSignalElectronLLDistributionGood->SetMaximum(1.35 * hSignalElectronLLDistributionGood->GetMaximum());
237
238 hSignalPionLLDistributionGood->GetXaxis()->SetTitle(PIDVarName + "PionLLSVDonly");
239 hSignalPionLLDistributionGood->GetYaxis()->SetTitle("Candidates, normalised");
240 hSignalPionLLDistributionGood->SetMaximum(1.35 * hSignalPionLLDistributionGood->GetMaximum());
241
242 hSignalKaonLLDistributionGood->GetXaxis()->SetTitle(PIDVarName + "KaonLLSVDonly");
243 hSignalKaonLLDistributionGood->GetYaxis()->SetTitle("Candidates, normalised");
244 hSignalKaonLLDistributionGood->SetMaximum(1.35 * hSignalKaonLLDistributionGood->GetMaximum());
245
246 hSignalProtonLLDistributionGood->GetXaxis()->SetTitle(PIDVarName + "ProtonLLSVDonly");
247 hSignalProtonLLDistributionGood->GetYaxis()->SetTitle("Candidates, normalised");
248 hSignalProtonLLDistributionGood->SetMaximum(1.35 * hSignalProtonLLDistributionGood->GetMaximum());
249
250 TCanvas* DistribCanvas = new TCanvas("DistribCanvas", "", 600, 600);
251 gPad->SetTopMargin(0.05);
252 gPad->SetRightMargin(0.05);
253 gPad->SetLeftMargin(0.13);
254 gPad->SetBottomMargin(0.12);
255
256 hSignalPIDDistribution->SetLineWidth(2);
257 hSignalPIDDistribution->SetLineColor(TColor::GetColor("#2166ac"));
258 hSignalPIDDistribution->Draw("hist");
259
260 DistribCanvas->Print("SVDdEdxValidation_Distribution_" + SignalVarNameFull + PIDVarName + PIDDetectorsName +
261 "_MomRange_" +
262 std::to_string(
263 MomLow)
264 .substr(0, 3) +
265 "_" + std::to_string(MomHigh).substr(0, 3) + ".pdf");
266
267 hSignalElectronLLDistribution->SetLineWidth(2);
268 hSignalPionLLDistribution->SetLineWidth(2);
269 hSignalKaonLLDistribution->SetLineWidth(2);
270 hSignalProtonLLDistribution->SetLineWidth(2);
271
272 hSignalElectronLLDistributionGood->SetLineWidth(2);
273 hSignalPionLLDistributionGood->SetLineWidth(2);
274 hSignalKaonLLDistributionGood->SetLineWidth(2);
275 hSignalProtonLLDistributionGood->SetLineWidth(2);
276
277 hSignalElectronLLDistributionGood->SetLineColor(kBlack);
278 hSignalPionLLDistributionGood->SetLineColor(kBlack);
279 hSignalKaonLLDistributionGood->SetLineColor(kBlack);
280 hSignalProtonLLDistributionGood->SetLineColor(kBlack);
281
282 hSignalElectronLLDistribution->SetTitle("ElectronLL (SVD), all tracks");
283 hSignalPionLLDistribution->SetTitle("PionLL (SVD), all tracks");
284 hSignalKaonLLDistribution->SetTitle("KaonLL (SVD), all tracks");
285 hSignalProtonLLDistribution->SetTitle("ProtonLL (SVD), all tracks");
286
287 hSignalElectronLLDistributionGood->SetTitle("ElectronLL (SVD), tracks with dEdx info");
288 hSignalPionLLDistributionGood->SetTitle("PionLL (SVD), tracks with dEdx info");
289 hSignalKaonLLDistributionGood->SetTitle("KaonLL (SVD), tracks with dEdx info");
290 hSignalProtonLLDistributionGood->SetTitle("ProtonLL (SVD), tracks with dEdx info");
291
292 hSignalElectronLLDistribution->GetXaxis()->SetTitleSize(0.04);
293 hSignalElectronLLDistribution->GetYaxis()->SetTitleSize(0.04);
294 hSignalElectronLLDistribution->GetXaxis()->SetTitleOffset(1.0);
295 hSignalElectronLLDistribution->GetYaxis()->SetTitleOffset(1.3);
296 hSignalElectronLLDistribution->GetYaxis()->SetLabelSize(0.04);
297 hSignalElectronLLDistribution->GetXaxis()->SetLabelSize(0.04);
298
299 hSignalPionLLDistribution->GetXaxis()->SetTitleSize(0.04);
300 hSignalPionLLDistribution->GetYaxis()->SetTitleSize(0.04);
301 hSignalPionLLDistribution->GetXaxis()->SetTitleOffset(1.0);
302 hSignalPionLLDistribution->GetYaxis()->SetTitleOffset(1.3);
303 hSignalPionLLDistribution->GetYaxis()->SetLabelSize(0.04);
304 hSignalPionLLDistribution->GetXaxis()->SetLabelSize(0.04);
305
306 hSignalKaonLLDistribution->GetXaxis()->SetTitleSize(0.04);
307 hSignalKaonLLDistribution->GetYaxis()->SetTitleSize(0.04);
308 hSignalKaonLLDistribution->GetXaxis()->SetTitleOffset(1.0);
309 hSignalKaonLLDistribution->GetYaxis()->SetTitleOffset(1.3);
310 hSignalKaonLLDistribution->GetYaxis()->SetLabelSize(0.04);
311 hSignalKaonLLDistribution->GetXaxis()->SetLabelSize(0.04);
312
313 hSignalProtonLLDistribution->GetXaxis()->SetTitleSize(0.04);
314 hSignalProtonLLDistribution->GetYaxis()->SetTitleSize(0.04);
315 hSignalProtonLLDistribution->GetXaxis()->SetTitleOffset(1.0);
316 hSignalProtonLLDistribution->GetYaxis()->SetTitleOffset(1.3);
317 hSignalProtonLLDistribution->GetYaxis()->SetLabelSize(0.04);
318 hSignalProtonLLDistribution->GetXaxis()->SetLabelSize(0.04);
319
320 TCanvas* LLCanvas = new TCanvas("LLCanvas", "", 900, 700);
321
322 gPad->SetTopMargin(0.05);
323 gPad->SetRightMargin(0.05);
324 gPad->SetLeftMargin(0.13);
325 gPad->SetBottomMargin(0.12);
326
327 LLCanvas->Divide(2, 2, 0.01, 0.01);
328 LLCanvas->cd(1);
329 hSignalElectronLLDistribution->Draw("hist");
330 LLCanvas->cd(2);
331 hSignalPionLLDistribution->Draw("hist");
332 LLCanvas->cd(3);
333 hSignalKaonLLDistribution->Draw("hist");
334 LLCanvas->cd(4);
335 hSignalProtonLLDistribution->Draw("hist");
336
337 TCanvas* LLCanvasGood = new TCanvas("LLCanvasGood", "", 900, 700);
338
339 gPad->SetTopMargin(0.05);
340 gPad->SetRightMargin(0.05);
341 gPad->SetLeftMargin(0.13);
342 gPad->SetBottomMargin(0.12);
343
344 LLCanvasGood->Divide(2, 2, 0.01, 0.01);
345 LLCanvasGood->cd(1);
346 hSignalElectronLLDistributionGood->Draw("hist");
347 LLCanvasGood->cd(2);
348 hSignalPionLLDistributionGood->Draw("hist");
349 LLCanvasGood->cd(3);
350 hSignalKaonLLDistributionGood->Draw("hist");
351 LLCanvasGood->cd(4);
352 hSignalProtonLLDistributionGood->Draw("hist");
353
354 LLCanvas->Print("SVDdEdxValidation_LLDistributions_" + SignalVarNameFull +
355 "_SVDonly_MomRange_" +
356 std::to_string(
357 MomLow)
358 .substr(0, 3) +
359 "_" + std::to_string(MomHigh).substr(0, 3) + ".pdf");
360
361 LLCanvasGood->Print("SVDdEdxValidation_LLDistributions_GoodSVDTracks_" + SignalVarNameFull +
362 "_SVDonly_MomRange_" +
363 std::to_string(
364 MomLow)
365 .substr(0, 3) +
366 "_" + std::to_string(MomHigh).substr(0, 3) + ".pdf");
367
368 TFile DistribFile("SVDdEdxValidation_Distribution_" + SignalVarNameFull + PIDVarName + PIDDetectorsName +
369 "_MomRange_" +
370 std::to_string(
371 MomLow)
372 .substr(0, 3) +
373 "_" + std::to_string(MomHigh).substr(0, 3) + ".root",
374 "RECREATE");
375 hSignalPIDDistribution->SetLineColor(kBlack);
376 hSignalPIDDistribution->Write();
377 DistribFile.Close();
378 delete DistribCanvas;
379
380 TFile LLDistribFile(TString("SVDdEdxValidation_LLDistributions_" + SignalVarNameFull + "_SVDonly_MomRange_" +
381 std::to_string(
382 MomLow)
383 .substr(0, 3) +
384 "_" + std::to_string(MomHigh).substr(0, 3) + ".root"),
385 "RECREATE");
386 hSignalElectronLLDistribution->Write();
387 hSignalPionLLDistribution->Write();
388 hSignalKaonLLDistribution->Write();
389 hSignalProtonLLDistribution->Write();
390 LLDistribFile.Close();
391 delete LLCanvas;
392 delete LLCanvasGood;
393 }
394
395 // ---------- Momentum distributions (for efficiency determination) ----------
396
397 SignalTree->Draw(Form("%sMomentum>>hAllSignal(%i,%f,%f)", SignalVarName.Data(), nbins, MomLow, MomHigh),
398 SignalWeightName + " * (" + SignalFiducialCut + ")", "goff");
399 SignalTree->Draw(Form("%sMomentum>>hSelectedSignal(%i,%f,%f)", SignalVarName.Data(), nbins, MomLow, MomHigh),
400 SignalWeightName + " * (" + SignalVarName + PIDVarName + PIDDetectorsName + ">" + PIDCut + "&&" + SignalFiducialCut +
401 ")",
402 "goff");
403
404 FakeTree->Draw(Form("%sMomentum>>hAllFakes(%i,%f,%f)", FakeVarName.Data(), nbins, MomLow, MomHigh),
405 FakeWeightName + " * (" + FakesFiducialCut + ")", "goff");
406 FakeTree->Draw(Form("%sMomentum>>hSelectedFakes(%i,%f,%f)", FakeVarName.Data(), nbins, MomLow, MomHigh),
407 FakeWeightName + " * (" + FakeVarName + PIDVarName + PIDDetectorsName + ">" + PIDCut + "&&" + FakesFiducialCut + ")",
408 "goff");
409
410 TH1F* hAllSignal = (TH1F*)gDirectory->Get("hAllSignal");
411 TH1F* hSelectedSignal = (TH1F*)gDirectory->Get("hSelectedSignal");
412 TH1F* hAllFakes = (TH1F*)gDirectory->Get("hAllFakes");
413 TH1F* hSelectedFakes = (TH1F*)gDirectory->Get("hSelectedFakes");
414
415 // ---------- Add slow pions to the pion dataset ----------
416 if (strncmp(SignalVarName.Data(), "PionD", 5) == 0) {
417 SignalTree->Draw(Form("SlowPionMomentum>>hAllSignalSlow(%i,%f,%f)", nbins, MomLow, MomHigh),
418 SignalWeightName + " * (" + SignalFiducialCut + ")", "goff");
419 SignalTree->Draw(Form("SlowPionMomentum>>hSelectedSignalSlow(%i,%f,%f)", nbins, MomLow, MomHigh),
420 SignalWeightName + " * (SlowPion" + PIDVarName + PIDDetectorsName + ">" + PIDCut + "&&" + SignalFiducialCut + ")", "goff");
421 TH1F* hAllSignalSlow = (TH1F*)gDirectory->Get("hAllSignalSlow");
422 TH1F* hSelectedSignalSlow = (TH1F*)gDirectory->Get("hSelectedSignalSlow");
423 hAllSignal->Add(hAllSignalSlow);
424 hSelectedSignal->Add(hSelectedSignalSlow);
425 }
426
427 if (strncmp(FakeVarName.Data(), "PionD", 5) == 0) {
428 FakeTree->Draw(Form("SlowPionMomentum>>hAllFakesSlow(%i,%f,%f)", nbins, MomLow, MomHigh),
429 FakeWeightName + " * (" + FakesFiducialCut + ")",
430 "goff");
431 FakeTree->Draw(Form("SlowPionMomentum>>hSelectedFakesSlow(%i,%f,%f)", nbins, MomLow, MomHigh),
432 FakeWeightName + " * (SlowPion" + PIDVarName + PIDDetectorsName + ">" + PIDCut + "&&" + FakesFiducialCut + ")", "goff");
433 TH1F* hAllFakesSlow = (TH1F*)gDirectory->Get("hAllFakesSlow");
434 TH1F* hSelectedFakesSlow = (TH1F*)gDirectory->Get("hSelectedFakesSlow");
435 hAllFakes->Add(hAllFakesSlow);
436 hSelectedFakes->Add(hSelectedFakesSlow);
437 }
438
439 TH1F* EffHistoSig = (TH1F*)hAllSignal->Clone("EffHistoSig"); // signal efficiency
440 TH1F* EffHistoFake = (TH1F*)hAllFakes->Clone("EffHistoFake"); // fakes efficiency
441
442 EffHistoSig->Divide(hSelectedSignal, hAllSignal);//, 1, 1, "B");
443 EffHistoFake->Divide(hSelectedFakes, hAllFakes);//, 1, 1, "B");
444
445 // PID plots
446 TH1F* hBase = new TH1F("hBase", "", 100, 0.0, MomHigh);
447 hBase->SetTitle(";Momentum [GeV];Efficiency");
448 hBase->SetMaximum(1.20);
449 hBase->SetMinimum(0.0);
450
451 TLegend* tleg1 = new TLegend(0.63, 0.82, 0.93, 0.94);
452 tleg1->AddEntry(EffHistoSig, SignalVarNameFull + " efficiency", "pl");
453 tleg1->AddEntry(EffHistoFake, FakeVarNameFull + " fake rate", "pl");
454
455 TCanvas* ResultCanvas = new TCanvas("ResultCanvas", "", 600, 600);
456 gPad->SetTopMargin(0.05);
457 gPad->SetRightMargin(0.05);
458 gPad->SetLeftMargin(0.13);
459 gPad->SetBottomMargin(0.12);
460
461 ResultCanvas->SetGrid();
462 hBase->Draw();
463 EffHistoSig->SetMarkerSize(1.5);
464 EffHistoSig->SetMarkerStyle(22);
465 EffHistoSig->SetMarkerColor(TColor::GetColor("#2166ac"));
466 EffHistoSig->SetLineColor(TColor::GetColor("#2166ac"));
467 EffHistoSig->Draw("P,same");
468
469 EffHistoFake->SetMarkerSize(1.5);
470 EffHistoFake->SetMarkerStyle(23);
471 EffHistoFake->SetMarkerColor(TColor::GetColor("#ef8a62"));
472 EffHistoFake->SetLineColor(TColor::GetColor("#ef8a62"));
473 EffHistoFake->Draw("P,same");
474
475 tleg1->Draw("same");
476
477 hBase->SetStats(0);
478 hBase->GetXaxis()->SetTitleSize(0.04);
479 hBase->GetYaxis()->SetTitleSize(0.04);
480 hBase->GetXaxis()->SetTitleOffset(1.0);
481 hBase->GetYaxis()->SetTitleOffset(1.3);
482 hBase->GetYaxis()->SetLabelSize(0.04);
483 hBase->GetXaxis()->SetLabelSize(0.04);
484
485 // std::setprecision(2);
486 ResultCanvas->Print("SVDdEdxValidation_Efficiency_" + SignalVarNameFull + "_vs_" + FakeVarNameFull + PIDVarName + "_" +
487 PIDDetectorsName +
488 "_Cut" +
489 PIDCut + "_MomRange_" + std::to_string(MomLow).substr(0, 3) + "_" + std::to_string(MomHigh).substr(0, 3) + ".pdf");
490 TFile ResultFile("SVDdEdxValidation_Efficiency_" + SignalVarNameFull + "_vs_" + FakeVarNameFull + PIDVarName + "_" +
491 PIDDetectorsName +
492 "_Cut" +
493 PIDCut + "_MomRange_" + std::to_string(MomLow).substr(0, 3) + "_" + std::to_string(MomHigh).substr(0, 3) + ".root",
494 "RECREATE");
495 EffHistoSig->SetLineColor(kBlack);
496 EffHistoSig->SetMarkerColor(kBlack);
497 EffHistoFake->SetLineColor(kBlack);
498 EffHistoFake->SetMarkerColor(kBlack);
499 EffHistoSig->Write();
500 EffHistoFake->Write();
501 ResultFile.Close();
502 delete ResultCanvas;
503 delete hBase;
504}
505
506void SVDdEdxValidationAlgorithm::PlotROCCurve(TTree* SignalTree, TString SignalWeightName, TString SignalVarName,
507 TString SignalVarNameFull, TTree* FakeTree, TString FakeWeightName, TString FakeVarName, TString FakeVarNameFull,
508 TString PIDVarName)
509{
510
511 if ((SignalTree == nullptr) || (FakeTree == nullptr)) {
512 B2FATAL("Invalid dataset, stopping here");
513 }
514
515 if ((SignalTree->GetEntries() == 0) || (FakeTree->GetEntries() == 0)) {
516 B2FATAL("The dataset is empty, stopping here");
517 }
518
519 if ((SignalTree->GetBranch(Form("%sMomentum", SignalVarName.Data())) == nullptr)
520 || (FakeTree->GetBranch(Form("%sMomentum", FakeVarName.Data())) == nullptr)) {
521 B2FATAL("Check the provided branch name, stopping here");
522 }
523
524 std::vector<TString> PIDDetectors;
525 PIDDetectors.clear();
526 PIDDetectors.push_back("ALL");
527 PIDDetectors.push_back("noSVD");
528
529 std::vector<double> SignalEfficiencyALL, FakeEfficiencyALL;
530 SignalEfficiencyALL.reserve(m_NumROCpoints);
531 FakeEfficiencyALL.reserve(m_NumROCpoints);
532 std::vector<double> SignalEfficiencynoSVD, FakeEfficiencynoSVD;
533 SignalEfficiencynoSVD.reserve(m_NumROCpoints);
534 FakeEfficiencynoSVD.reserve(m_NumROCpoints);
535
536 TString SignalFiducialCut = SignalVarName + PIDVarName + "noSVD>=0"; // sanity cuts to reject events with NaN
537 TString FakesFiducialCut = FakeVarName + PIDVarName + "noSVD>=0";
538 TString SignalFiducialCutSlow = "SlowPion" + PIDVarName + "noSVD>=0";
539 TString FakesFiducialCutSlow = "SlowPion" + PIDVarName + "noSVD>=0";
540
541 // calculate efficiencies
542
543 TCut AllSignalCut = SignalFiducialCut * Form("%sMomentum>%f && %sMomentum<%f", SignalVarName.Data(), m_MomLowROC,
544 SignalVarName.Data(), m_MomHighROC);
545
546 double AllSignalIntegral, SelectedSignalIntegral;
547
548 auto DataFrameSignalAll = RDataFrame(*SignalTree).Filter(AllSignalCut.GetTitle());
549
550 if (SignalWeightName == "1") {
551 AllSignalIntegral = DataFrameSignalAll.Count().GetValue();
552 } else {
553 AllSignalIntegral = DataFrameSignalAll.Sum(SignalWeightName).GetValue();
554 }
555
556 std::unique_ptr<ROOT::RDF::RNode> DataFrameSlowSignalAll;
557
558 if (strncmp(SignalVarName.Data(), "PionD", 5) == 0) {
559 TString SignalVarNameSlow = "SlowPion";
560 TCut AllSignalCutSlow = TCut(SignalFiducialCut) * TCut(SignalFiducialCutSlow) * Form("(%sMomentum>%f && %sMomentum<%f)",
561 SignalVarNameSlow.Data(), m_MomLowROC, SignalVarNameSlow.Data(), m_MomHighROC);
562 DataFrameSlowSignalAll = std::make_unique<ROOT::RDF::RNode>(RDataFrame(*SignalTree).Filter(AllSignalCutSlow.GetTitle()));
563
564 if (SignalWeightName == "1") {
565 AllSignalIntegral += DataFrameSlowSignalAll->Count().GetValue();
566 } else {
567 AllSignalIntegral += DataFrameSlowSignalAll->Sum(SignalWeightName).GetValue();
568 }
569 }
570
571 for (unsigned int i = 0; i < PIDDetectors.size(); i++) {
572 for (unsigned int j = 0; j < m_NumROCpoints; ++j) {
573 delete gROOT->FindObject("PIDCut");
574
575 // scan cut values from 0 to 1, with a denser scan closer to 0 or 1, to get a nicer ROC curve
576 double x = 1. / m_NumROCpoints * j;
577 TString PIDCut = TString::Format("%f", 1. / (1 + TMath::Power(x / (1 - x), -3)));
578
579 TCut SelectedSignalCut = Form("(%s%s%s > %s)", SignalVarName.Data(), PIDVarName.Data(), PIDDetectors[i].Data(), PIDCut.Data());
580
581 if (SignalWeightName == "1") {
582 SelectedSignalIntegral = DataFrameSignalAll.Filter(SelectedSignalCut.GetTitle()).Count().GetValue();
583 } else {
584 SelectedSignalIntegral = DataFrameSignalAll.Filter(SelectedSignalCut.GetTitle()).Sum(SignalWeightName).GetValue();
585 }
586
587 // special treatment for pions: add also the slow pions from Dstar to gain low-momentum coverage
588 if (strncmp(SignalVarName.Data(), "PionD", 5) == 0) {
589 TString SignalVarNameSlow = "SlowPion";
590 TCut SelectedSignalCutSlow = Form("(%s%s%s > %s)", SignalVarNameSlow.Data(), PIDVarName.Data(), PIDDetectors[i].Data(),
591 PIDCut.Data());
592
593 if (SignalWeightName == "1") {
594 SelectedSignalIntegral += DataFrameSlowSignalAll->Filter(SelectedSignalCutSlow.GetTitle()).Count().GetValue();
595 } else {
596 SelectedSignalIntegral += DataFrameSlowSignalAll->Filter(SelectedSignalCutSlow.GetTitle()).Sum(SignalWeightName).GetValue();
597 }
598 }
599
600 if (PIDDetectors[i] == "ALL") {
601 SignalEfficiencyALL.push_back(SelectedSignalIntegral / AllSignalIntegral);
602 }
603
604 if (PIDDetectors[i] == "noSVD") {
605 SignalEfficiencynoSVD.push_back(SelectedSignalIntegral / AllSignalIntegral);
606 }
607 }
608 }
609
610 // calculate fake rates
611
612 TCut AllFakeCut = FakesFiducialCut * Form("%sMomentum>%f && %sMomentum<%f", FakeVarName.Data(), m_MomLowROC, FakeVarName.Data(),
613 m_MomHighROC);
614
615 double AllFakeIntegral, SelectedFakeIntegral;
616 auto DataFrameFakeAll = RDataFrame(*FakeTree).Filter(AllFakeCut.GetTitle());
617
618 if (FakeWeightName == "1") {
619 AllFakeIntegral = DataFrameFakeAll.Count().GetValue();
620 } else {
621 AllFakeIntegral = DataFrameFakeAll.Sum(FakeWeightName).GetValue();
622 }
623
624 std::unique_ptr<ROOT::RDF::RNode> DataFrameSlowFakeAll;
625
626 // special treatment for pions: add also the slow pions from Dstar to gain low-momentum coverage
627 if (strncmp(FakeVarName.Data(), "PionD", 5) == 0) {
628
629 TString FakeVarNameSlow = "SlowPion";
630 TCut AllFakeCutSlow = TCut(FakesFiducialCut) * TCut(FakesFiducialCutSlow) * Form("(%sMomentum>%f && %sMomentum<%f)",
631 FakeVarNameSlow.Data(), m_MomLowROC, FakeVarNameSlow.Data(), m_MomHighROC);
632 DataFrameSlowFakeAll = std::make_unique<ROOT::RDF::RNode>(RDataFrame(*FakeTree).Filter(AllFakeCutSlow.GetTitle()));
633
634 if (FakeWeightName == "1") {
635 AllFakeIntegral += DataFrameSlowFakeAll->Count().GetValue();
636 } else {
637 AllFakeIntegral += DataFrameSlowFakeAll->Sum(FakeWeightName).GetValue();
638 }
639 }
640
641 for (unsigned int i = 0; i < PIDDetectors.size(); i++) {
642 for (unsigned int j = 0; j < m_NumROCpoints; ++j) {
643 delete gROOT->FindObject("PIDCut");
644 delete gROOT->FindObject("hAllFakes");
645 delete gROOT->FindObject("hSelectedFakes");
646
647 // scan cut values from 0 to 1, with a denser scan closer to 0 or 1, to get a nicer ROC curve
648 double x = 1. / m_NumROCpoints * j;
649 TString PIDCut = TString::Format("%f", 1. / (1 + TMath::Power(x / (1 - x), -3)));
650
651 TCut SelectedFakeCut = Form("(%s%s%s > %s)", FakeVarName.Data(), PIDVarName.Data(), PIDDetectors[i].Data(), PIDCut.Data());
652
653 if (FakeWeightName == "1") {
654 SelectedFakeIntegral = DataFrameFakeAll.Filter(SelectedFakeCut.GetTitle()).Count().GetValue();
655 } else {
656 SelectedFakeIntegral = DataFrameFakeAll.Filter(SelectedFakeCut.GetTitle()).Sum(FakeWeightName).GetValue();
657 }
658
659 if (strncmp(FakeVarName.Data(), "PionD", 5) == 0) {
660 TString FakeVarNameSlow = "SlowPion";
661
662 TCut SelectedFakeCutSlow = Form("(%s%s%s > %s)", FakeVarNameSlow.Data(), PIDVarName.Data(), PIDDetectors[i].Data(), PIDCut.Data());
663
664 if (FakeWeightName == "1") {
665 SelectedFakeIntegral += DataFrameSlowFakeAll->Filter(SelectedFakeCutSlow.GetTitle()).Count().GetValue();
666 } else {
667 SelectedFakeIntegral += DataFrameSlowFakeAll->Filter(SelectedFakeCutSlow.GetTitle()).Sum(FakeWeightName).GetValue();
668 }
669 }
670
671 if (PIDDetectors[i] == "ALL") {
672 FakeEfficiencyALL.push_back(SelectedFakeIntegral / AllFakeIntegral);
673 }
674
675 if (PIDDetectors[i] == "noSVD") {
676 FakeEfficiencynoSVD.push_back(SelectedFakeIntegral / AllFakeIntegral);
677 }
678 }
679 }
680
681 auto ResultCanvas = new TCanvas("ResultCanvas", "", 600, 400);
682 TMultiGraph* hmgraph = new TMultiGraph();
683
684 // efficiency vs fake rate graph
685 TGraph* hgraphALL = new TGraph(m_NumROCpoints, FakeEfficiencyALL.data(), SignalEfficiencyALL.data());
686 hgraphALL->SetMarkerColor(TColor::GetColor("#2166ac"));
687 hgraphALL->SetMarkerStyle(20);
688 hgraphALL->SetLineColor(TColor::GetColor("#2166ac"));
689 hgraphALL->SetLineWidth(3);
690 hgraphALL->SetDrawOption("AP*");
691 hgraphALL->SetTitle("with SVD");
692
693 TGraph* hgraphnoSVD = new TGraph(m_NumROCpoints, FakeEfficiencynoSVD.data(), SignalEfficiencynoSVD.data());
694 hgraphnoSVD->SetMarkerColor(TColor::GetColor("#ef8a62"));
695 hgraphnoSVD->SetLineColor(TColor::GetColor("#ef8a62"));
696 hgraphnoSVD->SetLineWidth(3);
697 hgraphnoSVD->SetMarkerStyle(22);
698 hgraphnoSVD->SetDrawOption("P*");
699 hgraphnoSVD->SetTitle("without SVD");
700
701 hmgraph->Add(hgraphALL);
702 hmgraph->Add(hgraphnoSVD);
703 hmgraph->Draw("A");
704 hmgraph->GetHistogram()->GetXaxis()->SetTitle(FakeVarNameFull + " fake rate");
705 hmgraph->GetHistogram()->GetYaxis()->SetTitle(SignalVarNameFull + " signal efficiency");
706
707 ResultCanvas->BuildLegend(0.6, 0.25, 0.9, 0.5);
708 ResultCanvas->SetGrid();
709
710 ResultCanvas->Print("SVDdEdxValidation_ROC_curve_" + SignalVarNameFull + "_vs_" + FakeVarNameFull + PIDVarName + "_MomRange" +
711 std::to_string(m_MomLowROC).substr(0, 3) + "_" + std::to_string(m_MomHighROC).substr(0, 3) + ".pdf");
712
713 TFile ResultFile("SVDdEdxValidation_ROC_curve_" + SignalVarNameFull + "_vs_" + FakeVarNameFull + PIDVarName + "_MomRange" +
714 std::to_string(m_MomLowROC).substr(0, 3) + "_" + std::to_string(m_MomHighROC).substr(0, 3) + ".root",
715 "RECREATE");
716 hmgraph->Write();
717 ResultFile.Close();
718
719 delete ResultCanvas;
720}
721
722TTree* SVDdEdxValidationAlgorithm::LambdaMassFit(std::shared_ptr<TTree> preselTree)
723{
724 B2INFO("Configuring the Lambda fit...");
725 gROOT->SetBatch(true);
726 RooMsgService::instance().setGlobalKillBelow(RooFit::WARNING);
727
728 RooRealVar InvM("InvM", "m(p^{+}#pi^{-})", 1.1, 1.13, "GeV/c^{2}");
729
730 RooRealVar ProtonMomentum("ProtonMomentum", "momentum for p", -1.e8, 1.e8);
731 RooRealVar ProtonSVDdEdx("ProtonSVDdEdx", "", -1.e8, 1.e8);
732
733 RooRealVar exp("exp", "experiment number", 0, 1.e5);
734 RooRealVar run("run", "run number", 0, 1.e7);
735
736 RooRealVar ProtonProtonIDALL("ProtonProtonIDALL", "", -1.e8, 1.e8);
737 RooRealVar ProtonProtonIDSVDonly("ProtonProtonIDSVDonly", "", -1.e8, 1.e8);
738 RooRealVar ProtonProtonIDnoSVD("ProtonProtonIDnoSVD", "", -1.e8, 1.e8);
739
740 RooRealVar ProtonElectronLLSVDonly("ProtonElectronLLSVDonly", "", -1.e8, 1.e8);
741 RooRealVar ProtonPionLLSVDonly("ProtonPionLLSVDonly", "", -1.e8, 1.e8);
742 RooRealVar ProtonKaonLLSVDonly("ProtonKaonLLSVDonly", "", -1.e8, 1.e8);
743 RooRealVar ProtonProtonLLSVDonly("ProtonProtonLLSVDonly", "", -1.e8, 1.e8);
744
745 RooRealVar ProtonBinaryProtonPionIDALL("ProtonBinaryProtonPionIDALL", "", -1.e8, 1.e8);
746 RooRealVar ProtonBinaryProtonPionIDSVDonly("ProtonBinaryProtonPionIDSVDonly", "", -1.e8, 1.e8);
747 RooRealVar ProtonBinaryProtonPionIDnoSVD("ProtonBinaryProtonPionIDnoSVD", "", -1.e8, 1.e8);
748
749 RooRealVar ProtonBinaryProtonKaonIDALL("ProtonBinaryProtonKaonIDALL", "", -1.e8, 1.e8);
750 RooRealVar ProtonBinaryProtonKaonIDSVDonly("ProtonBinaryProtonKaonIDSVDonly", "", -1.e8, 1.e8);
751 RooRealVar ProtonBinaryProtonKaonIDnoSVD("ProtonBinaryProtonKaonIDnoSVD", "", -1.e8, 1.e8);
752
753 RooRealVar ProtonBinaryProtonElectronIDALL("ProtonBinaryProtonElectronIDALL", "", -1.e8, 1.e8);
754 RooRealVar ProtonBinaryProtonElectronIDSVDonly("ProtonBinaryProtonElectronIDSVDonly", "", -1.e8, 1.e8);
755 RooRealVar ProtonBinaryProtonElectronIDnoSVD("ProtonBinaryProtonElectronIDnoSVD", "", -1.e8, 1.e8);
756
757 RooRealVar ProtonBinaryPionProtonIDALL("ProtonBinaryPionProtonIDALL", "", -1.e8, 1.e8);
758 RooRealVar ProtonBinaryPionProtonIDSVDonly("ProtonBinaryPionProtonIDSVDonly", "", -1.e8, 1.e8);
759 RooRealVar ProtonBinaryPionProtonIDnoSVD("ProtonBinaryPionProtonIDnoSVD", "", -1.e8, 1.e8);
760
761 RooRealVar ProtonBinaryKaonProtonIDALL("ProtonBinaryKaonProtonIDALL", "", -1.e8, 1.e8);
762 RooRealVar ProtonBinaryKaonProtonIDSVDonly("ProtonBinaryKaonProtonIDSVDonly", "", -1.e8, 1.e8);
763 RooRealVar ProtonBinaryKaonProtonIDnoSVD("ProtonBinaryKaonProtonIDnoSVD", "", -1.e8, 1.e8);
764
765 RooRealVar ProtonBinaryElectronProtonIDALL("ProtonBinaryElectronProtonIDALL", "", -1.e8, 1.e8);
766 RooRealVar ProtonBinaryElectronProtonIDSVDonly("ProtonBinaryElectronProtonIDSVDonly", "", -1.e8, 1.e8);
767 RooRealVar ProtonBinaryElectronProtonIDnoSVD("ProtonBinaryElectronProtonIDnoSVD", "", -1.e8, 1.e8);
768
769 auto variables = new RooArgSet();
770
771 variables->add(InvM);
772
773 variables->add(ProtonMomentum);
774 variables->add(ProtonSVDdEdx);
775 variables->add(exp);
776 variables->add(run);
777
778 variables->add(ProtonProtonIDALL);
779 variables->add(ProtonProtonIDSVDonly);
780 variables->add(ProtonProtonIDnoSVD);
781 variables->add(ProtonElectronLLSVDonly);
782 variables->add(ProtonPionLLSVDonly);
783 variables->add(ProtonKaonLLSVDonly);
784 variables->add(ProtonProtonLLSVDonly);
785 variables->add(ProtonBinaryProtonPionIDALL);
786 variables->add(ProtonBinaryProtonPionIDSVDonly);
787 variables->add(ProtonBinaryProtonPionIDnoSVD);
788 variables->add(ProtonBinaryProtonKaonIDALL);
789 variables->add(ProtonBinaryProtonKaonIDSVDonly);
790 variables->add(ProtonBinaryProtonKaonIDnoSVD);
791 variables->add(ProtonBinaryProtonElectronIDALL);
792 variables->add(ProtonBinaryProtonElectronIDSVDonly);
793 variables->add(ProtonBinaryProtonElectronIDnoSVD);
794 variables->add(ProtonBinaryPionProtonIDALL);
795 variables->add(ProtonBinaryPionProtonIDSVDonly);
796 variables->add(ProtonBinaryPionProtonIDnoSVD);
797 variables->add(ProtonBinaryKaonProtonIDALL);
798 variables->add(ProtonBinaryKaonProtonIDSVDonly);
799 variables->add(ProtonBinaryKaonProtonIDnoSVD);
800 variables->add(ProtonBinaryElectronProtonIDALL);
801 variables->add(ProtonBinaryElectronProtonIDSVDonly);
802 variables->add(ProtonBinaryElectronProtonIDnoSVD);
803
804 RooDataSet* LambdaDataset = new RooDataSet("LambdaDataset", "LambdaDataset", preselTree.get(), *variables);
805
806 if (LambdaDataset->sumEntries() == 0) {
807 B2FATAL("The Lambda dataset is empty, stopping here");
808 }
809
810 // the signal PDF; might be revisited at a later point
811
812 RooRealVar GaussMean("GaussMean", " GaussMean", 1.116, 1.111, 1.12);
813 RooRealVar GaussSigma("GaussSigma", "#sigma_{1}", 3.e-3, 3.e-5, 10.e-3);
814 RooGaussian LambdaGauss("LambdaGauss", "LambdaGauss", InvM, GaussMean, GaussSigma);
815
816 /* temporary RooRealVar sigmaBifurGaussL1 and sigmaBifurGaussR1 to replace
817 * RooRealVar resolutionParamL("resolutionParamL", "resolutionParamL", 0.4, 5.e-4, 1.0);
818 * RooRealVar resolutionParamR("resolutionParamR", "resolutionParamR", 0.4, 5.e-4, 1.0);
819 * RooFormulaVar sigmaBifurGaussL1("sigmaBifurGaussL1", "resolutionParamL*GaussSigma", RooArgSet(resolutionParamL, GaussSigma));
820 * RooFormulaVar sigmaBifurGaussR1("sigmaBifurGaussR1", "resolutionParamR*GaussSigma", RooArgSet(resolutionParamR, GaussSigma));
821 */
822 RooRealVar sigmaBifurGaussL1("sigmaBifurGaussL1", "sigma left", 0.4 * 3.e-3, 3.e-5, 10.e-3);
823 RooRealVar sigmaBifurGaussR1("sigmaBifurGaussR1", "sigma right", 0.4 * 3.e-3, 3.e-5, 10.e-3);
824 RooBifurGauss LambdaBifurGauss("LambdaBifurGauss", "LambdaBifurGauss", InvM, GaussMean, sigmaBifurGaussL1, sigmaBifurGaussR1);
825
826 /* temporary RooRealVar sigmaBifurGaussL2 to replace
827 * RooRealVar resolutionParam2("resolutionParam2", "resolutionParam2", 0.2, 5.e-4, 1.0);
828 * sigmaBifurGaussL2("sigmaBifurGaussL2", "resolutionParam2*GaussSigma", RooArgSet(resolutionParam2, GaussSigma));
829 */
830 RooRealVar sigmaBifurGaussL2("sigmaBifurGaussL2", "sigmaBifurGaussL2", 0.2 * 3.e-3, 3.e-5, 10.e-3);
831 RooGaussian LambdaBifurGauss2("LambdaBifurGauss2", "LambdaBifurGauss2", InvM, GaussMean, sigmaBifurGaussL2);
832
833 RooRealVar fracBifurGaussYield("fracBifurGaussYield", "fracBifurGaussYield", 0.3, 5.e-4, 1.0);
834 RooRealVar fracGaussYield("fracGaussYield", "fracGaussYield", 0.8, 5.e-4, 1.0);
835
836 RooAddPdf LambdaCombinedBifurGauss("LambdaCombinedBifurGauss", "LambdaBifurGauss + LambdaBifurGauss2 ", RooArgList(LambdaBifurGauss,
837 LambdaBifurGauss2), RooArgList(fracBifurGaussYield));
838
839 RooAddPdf LambdaSignalPDF("LambdaSignalPDF", "LambdaCombinedBifurGauss + LambdaGauss", RooArgList(LambdaCombinedBifurGauss,
840 LambdaGauss), RooArgList(fracGaussYield));
841
842 // Background PDF
843 RooRealVar BkgPolyCoef0("BkgPolyCoef0", "BkgPolyCoef0", 0.1, 0., 1.5);
844 RooRealVar BkgPolyCoef1("BkgPolyCoef1", "BkgPolyCoef1", -0.5, -1.5, -1.e-3);
845 RooChebychev BkgPolyPDF("BkgPolyPDF", "BkgPolyPDF", InvM, RooArgList(BkgPolyCoef0, BkgPolyCoef1));
846
847 RooRealVar nSignalLambda("nSignalLambda", "nSignalLambda", 0.6 * preselTree->GetEntries(), 0., 0.99 * preselTree->GetEntries());
848 RooRealVar nBkgLambda("nBkgLambda", "nBkgLambda", 0.4 * preselTree->GetEntries(), 0., 0.99 * preselTree->GetEntries());
849 RooAddPdf totalPDFLambda("totalPDFLambda", "totalPDFLambda pdf", RooArgList(LambdaSignalPDF, BkgPolyPDF),
850 RooArgList(nSignalLambda, nBkgLambda));
851
852 B2INFO("Lambda: Start fitting...");
853 RooFitResult* LambdaFitResult = totalPDFLambda.fitTo(*LambdaDataset, Save(kTRUE), PrintLevel(-1));
854
855 int status = LambdaFitResult->status();
856 int covqual = LambdaFitResult->covQual();
857 double diff = nSignalLambda.getValV() + nBkgLambda.getValV() - LambdaDataset->sumEntries();
858
859 B2INFO("Lambda: Fit status: " << status << "; covariance quality: " << covqual);
860 // if the fit is not healthy, try again once before giving up, with a slightly different setup:
861 if ((status > 0) || (TMath::Abs(diff) > 1.) || (nSignalLambda.getError() < sqrt(nSignalLambda.getValV()))
862 || (nSignalLambda.getError() > (nSignalLambda.getValV()))) {
863
864 LambdaFitResult = totalPDFLambda.fitTo(*LambdaDataset, Save(), Strategy(2), Offset(1));
865 status = LambdaFitResult->status();
866 covqual = LambdaFitResult->covQual();
867 diff = nSignalLambda.getValV() + nBkgLambda.getValV() - LambdaDataset->sumEntries();
868 }
869
870 if ((status > 0) || (TMath::Abs(diff) > 1.) || (nSignalLambda.getError() < sqrt(nSignalLambda.getValV()))
871 || (nSignalLambda.getError() > (nSignalLambda.getValV()))) {
872 B2WARNING("Lambda: Fit problem: fit status " << status << "; sum of component yields minus the dataset yield is " << diff <<
873 "; signal yield is " << nSignalLambda.getValV() << ", while its uncertainty is " << nSignalLambda.getError());
874 }
875 if (covqual < 2) {
876 B2INFO("Lambda: Fit warning: covariance quality " << covqual);
877 }
878
879 TCanvas* canvLambda = new TCanvas("canvLambda", "canvLambda");
880 RooPlot* LambdaFitFrame = LambdaDataset->plotOn(InvM.frame(130));
881 totalPDFLambda.plotOn(LambdaFitFrame, LineColor(TColor::GetColor("#4575b4")));
882
883 double chisquare = LambdaFitFrame->chiSquare();
884 B2INFO("Lambda: Fit chi2 = " << chisquare);
885 totalPDFLambda.paramOn(LambdaFitFrame, Layout(0.6, 0.96, 0.93), Format("NEU", AutoPrecision(2)));
886 LambdaFitFrame->getAttText()->SetTextSize(0.03);
887
888 totalPDFLambda.plotOn(LambdaFitFrame, Components("LambdaSignalPDF"), LineColor(TColor::GetColor("#d73027")));
889 totalPDFLambda.plotOn(LambdaFitFrame, Components("BkgPolyPDF"), LineColor(TColor::GetColor("#fc8d59")));
890 totalPDFLambda.plotOn(LambdaFitFrame, LineColor(TColor::GetColor("#4575b4")));
891
892 LambdaFitFrame->GetXaxis()->SetTitle("m(p#pi^{-}) (GeV/c^{2})");
893
894 LambdaFitFrame->Draw();
895
896 if (m_isMakePlots) {
897 canvLambda->Print("SVDdEdxValidationFitLambda.pdf");
898 TFile LambdaFitPlotFile("SVDdEdxValidationLambdaFitPlotFile.root", "RECREATE");
899 canvLambda->Write();
900 LambdaFitPlotFile.Close();
901 }
902 RooStats::SPlot* sPlotDatasetLambda = new RooStats::SPlot("sData", "An SPlot", *LambdaDataset, &totalPDFLambda,
903 RooArgList(nSignalLambda, nBkgLambda));
904
905 for (int iEvt = 0; iEvt < 5; iEvt++) {
906 if (TMath::Abs(sPlotDatasetLambda->GetSWeight(iEvt, "nSignalLambda") + sPlotDatasetLambda->GetSWeight(iEvt,
907 "nBkgLambda") - 1) > 5.e-3)
908 B2FATAL("Lambda: sPlot error: sum of weights not equal to 1");
909 }
910
911 RooDataSet* LambdaDatasetSWeighted = new RooDataSet(LambdaDataset->GetName(), LambdaDataset->GetTitle(), LambdaDataset,
912 *LambdaDataset->get());
913
914 TTree* treeLambda_sw = LambdaDatasetSWeighted->GetClonedTree();
915 treeLambda_sw->SetName("treeLambda_sw");
916
917 B2INFO("Lambda: sPlot done. ");
918
919 return treeLambda_sw;
920}
921
922TTree* SVDdEdxValidationAlgorithm::DstarMassFit(std::shared_ptr<TTree> preselTree)
923{
924 B2INFO("Configuring the Dstar fit...");
925 gROOT->SetBatch(true);
926 RooMsgService::instance().setGlobalKillBelow(RooFit::WARNING);
927
928 RooRealVar deltaM("deltaM", "m(D*)-m(D^{0})", 0.139545, 0.151, "GeV/c^{2}");
929
930 RooRealVar KaonMomentum("KaonMomentum", "momentum for Kaon(GeV)", -1.e8, 1.e8);
931 RooRealVar KaonSVDdEdx("KaonSVDdEdx", "", -1.e8, 1.e8);
932 RooRealVar PionDMomentum("PionDMomentum", "momentum for pion(GeV)", -1.e8, 1.e8);
933 RooRealVar PionDSVDdEdx("PionDSVDdEdx", "", -1.e8, 1.e8);
934 RooRealVar SlowPionMomentum("SlowPionMomentum", "momentum for slow pion(GeV)", -1.e8, 1.e8);
935 RooRealVar SlowPionSVDdEdx("SlowPionSVDdEdx", "", -1.e8, 1.e8);
936
937 RooRealVar exp("exp", "experiment number", 0, 1.e5);
938 RooRealVar run("run", "run number", 0, 1.e8);
939
940 RooRealVar KaonKaonIDALL("KaonKaonIDALL", "", -1.e8, 1.e8);
941 RooRealVar KaonKaonIDSVDonly("KaonKaonIDSVDonly", "", -1.e8, 1.e8);
942 RooRealVar KaonKaonIDnoSVD("KaonKaonIDnoSVD", "", -1.e8, 1.e8);
943
944 RooRealVar KaonPionIDALL("KaonPionIDALL", "", -1.e8, 1.e8);
945 RooRealVar KaonPionIDSVDonly("KaonPionIDSVDonly", "", -1.e8, 1.e8);
946 RooRealVar KaonPionIDnoSVD("KaonPionIDnoSVD", "", -1.e8, 1.e8);
947
948 RooRealVar KaonProtonIDALL("KaonProtonIDALL", "", -1.e8, 1.e8);
949 RooRealVar KaonProtonIDSVDonly("KaonProtonIDSVDonly", "", -1.e8, 1.e8);
950 RooRealVar KaonProtonIDnoSVD("KaonProtonIDnoSVD", "", -1.e8, 1.e8);
951
952 RooRealVar KaonElectronIDALL("KaonElectronIDALL", "", -1.e8, 1.e8);
953 RooRealVar KaonElectronIDSVDonly("KaonElectronIDSVDonly", "", -1.e8, 1.e8);
954 RooRealVar KaonElectronIDnoSVD("KaonElectronIDnoSVD", "", -1.e8, 1.e8);
955
956 RooRealVar KaonElectronLLSVDonly("KaonElectronLLSVDonly", "", -1.e8, 1.e8);
957 RooRealVar KaonPionLLSVDonly("KaonPionLLSVDonly", "", -1.e8, 1.e8);
958 RooRealVar KaonKaonLLSVDonly("KaonKaonLLSVDonly", "", -1.e8, 1.e8);
959 RooRealVar KaonProtonLLSVDonly("KaonProtonLLSVDonly", "", -1.e8, 1.e8);
960
961 RooRealVar KaonBinaryKaonPionIDALL("KaonBinaryKaonPionIDALL", "", -1.e8, 1.e8);
962 RooRealVar KaonBinaryKaonPionIDSVDonly("KaonBinaryKaonPionIDSVDonly", "", -1.e8, 1.e8);
963 RooRealVar KaonBinaryKaonPionIDnoSVD("KaonBinaryKaonPionIDnoSVD", "", -1.e8, 1.e8);
964
965 RooRealVar KaonBinaryPionKaonIDALL("KaonBinaryPionKaonIDALL", "", -1.e8, 1.e8);
966 RooRealVar KaonBinaryPionKaonIDSVDonly("KaonBinaryPionKaonIDSVDonly", "", -1.e8, 1.e8);
967 RooRealVar KaonBinaryPionKaonIDnoSVD("KaonBinaryPionKaonIDnoSVD", "", -1.e8, 1.e8);
968
969 RooRealVar KaonBinaryProtonKaonIDALL("KaonBinaryProtonKaonIDALL", "", -1.e8, 1.e8);
970 RooRealVar KaonBinaryProtonKaonIDSVDonly("KaonBinaryProtonKaonIDSVDonly", "", -1.e8, 1.e8);
971 RooRealVar KaonBinaryProtonKaonIDnoSVD("KaonBinaryProtonKaonIDnoSVD", "", -1.e8, 1.e8);
972
973 RooRealVar KaonBinaryElectronKaonIDALL("KaonBinaryElectronKaonIDALL", "", -1.e8, 1.e8);
974 RooRealVar KaonBinaryElectronKaonIDSVDonly("KaonBinaryElectronKaonIDSVDonly", "", -1.e8, 1.e8);
975 RooRealVar KaonBinaryElectronKaonIDnoSVD("KaonBinaryElectronKaonIDnoSVD", "", -1.e8, 1.e8);
976
977 RooRealVar PionDKaonIDALL("PionDKaonIDALL", "", -1.e8, 1.e8);
978 RooRealVar PionDKaonIDSVDonly("PionDKaonIDSVDonly", "", -1.e8, 1.e8);
979 RooRealVar PionDKaonIDnoSVD("PionDKaonIDnoSVD", "", -1.e8, 1.e8);
980
981 RooRealVar PionDPionIDALL("PionDPionIDALL", "", -1.e8, 1.e8);
982 RooRealVar PionDPionIDSVDonly("PionDPionIDSVDonly", "", -1.e8, 1.e8);
983 RooRealVar PionDPionIDnoSVD("PionDPionIDnoSVD", "", -1.e8, 1.e8);
984
985 RooRealVar PionDElectronIDALL("PionDElectronIDALL", "", -1.e8, 1.e8);
986 RooRealVar PionDElectronIDSVDonly("PionDElectronIDSVDonly", "", -1.e8, 1.e8);
987 RooRealVar PionDElectronIDnoSVD("PionDElectronIDnoSVD", "", -1.e8, 1.e8);
988
989 RooRealVar PionDProtonIDALL("PionDProtonIDALL", "", -1.e8, 1.e8);
990 RooRealVar PionDProtonIDSVDonly("PionDProtonIDSVDonly", "", -1.e8, 1.e8);
991 RooRealVar PionDProtonIDnoSVD("PionDProtonIDnoSVD", "", -1.e8, 1.e8);
992
993 RooRealVar PionDElectronLLSVDonly("PionDElectronLLSVDonly", "", -1.e8, 1.e8);
994 RooRealVar PionDPionLLSVDonly("PionDPionLLSVDonly", "", -1.e8, 1.e8);
995 RooRealVar PionDKaonLLSVDonly("PionDKaonLLSVDonly", "", -1.e8, 1.e8);
996 RooRealVar PionDProtonLLSVDonly("PionDProtonLLSVDonly", "", -1.e8, 1.e8);
997
998 RooRealVar PionDBinaryPionKaonIDALL("PionDBinaryPionKaonIDALL", "", -1.e8, 1.e8);
999 RooRealVar PionDBinaryPionKaonIDSVDonly("PionDBinaryPionKaonIDSVDonly", "", -1.e8, 1.e8);
1000 RooRealVar PionDBinaryPionKaonIDnoSVD("PionDBinaryPionKaonIDnoSVD", "", -1.e8, 1.e8);
1001
1002 RooRealVar PionDBinaryKaonPionIDALL("PionDBinaryKaonPionIDALL", "", -1.e8, 1.e8);
1003 RooRealVar PionDBinaryKaonPionIDSVDonly("PionDBinaryKaonPionIDSVDonly", "", -1.e8, 1.e8);
1004 RooRealVar PionDBinaryKaonPionIDnoSVD("PionDBinaryKaonPionIDnoSVD", "", -1.e8, 1.e8);
1005
1006 RooRealVar PionDBinaryProtonPionIDALL("PionDBinaryProtonPionIDALL", "", -1.e8, 1.e8);
1007 RooRealVar PionDBinaryProtonPionIDSVDonly("PionDBinaryProtonPionIDSVDonly", "", -1.e8, 1.e8);
1008 RooRealVar PionDBinaryProtonPionIDnoSVD("PionDBinaryProtonPionIDnoSVD", "", -1.e8, 1.e8);
1009
1010 RooRealVar PionDBinaryElectronPionIDALL("PionDBinaryElectronPionIDALL", "", -1.e8, 1.e8);
1011 RooRealVar PionDBinaryElectronPionIDSVDonly("PionDBinaryElectronPionIDSVDonly", "", -1.e8, 1.e8);
1012 RooRealVar PionDBinaryElectronPionIDnoSVD("PionDBinaryElectronPionIDnoSVD", "", -1.e8, 1.e8);
1013
1014 RooRealVar SlowPionKaonIDALL("SlowPionKaonIDALL", "", -1.e8, 1.e8);
1015 RooRealVar SlowPionKaonIDSVDonly("SlowPionKaonIDSVDonly", "", -1.e8, 1.e8);
1016 RooRealVar SlowPionKaonIDnoSVD("SlowPionKaonIDnoSVD", "", -1.e8, 1.e8);
1017
1018 RooRealVar SlowPionPionIDALL("SlowPionPionIDALL", "", -1.e8, 1.e8);
1019 RooRealVar SlowPionPionIDSVDonly("SlowPionPionIDSVDonly", "", -1.e8, 1.e8);
1020 RooRealVar SlowPionPionIDnoSVD("SlowPionPionIDnoSVD", "", -1.e8, 1.e8);
1021
1022 RooRealVar SlowPionElectronIDALL("SlowPionElectronIDALL", "", -1.e8, 1.e8);
1023 RooRealVar SlowPionElectronIDSVDonly("SlowPionElectronIDSVDonly", "", -1.e8, 1.e8);
1024 RooRealVar SlowPionElectronIDnoSVD("SlowPionElectronIDnoSVD", "", -1.e8, 1.e8);
1025
1026 RooRealVar SlowPionProtonIDALL("SlowPionProtonIDALL", "", -1.e8, 1.e8);
1027 RooRealVar SlowPionProtonIDSVDonly("SlowPionProtonIDSVDonly", "", -1.e8, 1.e8);
1028 RooRealVar SlowPionProtonIDnoSVD("SlowPionProtonIDnoSVD", "", -1.e8, 1.e8);
1029
1030 RooRealVar SlowPionElectronLLSVDonly("SlowPionElectronLLSVDonly", "", -1.e8, 1.e8);
1031 RooRealVar SlowPionPionLLSVDonly("SlowPionPionLLSVDonly", "", -1.e8, 1.e8);
1032 RooRealVar SlowPionKaonLLSVDonly("SlowPionKaonLLSVDonly", "", -1.e8, 1.e8);
1033 RooRealVar SlowPionProtonLLSVDonly("SlowPionProtonLLSVDonly", "", -1.e8, 1.e8);
1034
1035 RooRealVar SlowPionBinaryPionKaonIDALL("SlowPionBinaryPionKaonIDALL", "", -1.e8, 1.e8);
1036 RooRealVar SlowPionBinaryPionKaonIDSVDonly("SlowPionBinaryPionKaonIDSVDonly", "", -1.e8, 1.e8);
1037 RooRealVar SlowPionBinaryPionKaonIDnoSVD("SlowPionBinaryPionKaonIDnoSVD", "", -1.e8, 1.e8);
1038
1039 RooRealVar SlowPionBinaryKaonPionIDALL("SlowPionBinaryKaonPionIDALL", "", -1.e8, 1.e8);
1040 RooRealVar SlowPionBinaryKaonPionIDSVDonly("SlowPionBinaryKaonPionIDSVDonly", "", -1.e8, 1.e8);
1041 RooRealVar SlowPionBinaryKaonPionIDnoSVD("SlowPionBinaryKaonPionIDnoSVD", "", -1.e8, 1.e8);
1042
1043 RooRealVar SlowPionBinaryProtonPionIDALL("SlowPionBinaryProtonPionIDALL", "", -1.e8, 1.e8);
1044 RooRealVar SlowPionBinaryProtonPionIDSVDonly("SlowPionBinaryProtonPionIDSVDonly", "", -1.e8, 1.e8);
1045 RooRealVar SlowPionBinaryProtonPionIDnoSVD("SlowPionBinaryProtonPionIDnoSVD", "", -1.e8, 1.e8);
1046
1047 RooRealVar SlowPionBinaryElectronPionIDALL("SlowPionBinaryElectronPionIDALL", "", -1.e8, 1.e8);
1048 RooRealVar SlowPionBinaryElectronPionIDSVDonly("SlowPionBinaryElectronPionIDSVDonly", "", -1.e8, 1.e8);
1049 RooRealVar SlowPionBinaryElectronPionIDnoSVD("SlowPionBinaryElectronPionIDnoSVD", "", -1.e8, 1.e8);
1050
1051 auto variables = new RooArgSet();
1052 variables->add(deltaM);
1053 variables->add(KaonMomentum);
1054 variables->add(KaonSVDdEdx);
1055 variables->add(PionDMomentum);
1056 variables->add(PionDSVDdEdx);
1057 variables->add(SlowPionMomentum);
1058 variables->add(SlowPionSVDdEdx);
1059 variables->add(exp);
1060 variables->add(run);
1061
1062 variables->add(KaonKaonIDALL);
1063 variables->add(KaonKaonIDSVDonly);
1064 variables->add(KaonKaonIDnoSVD);
1065 variables->add(KaonPionIDALL);
1066 variables->add(KaonPionIDSVDonly);
1067 variables->add(KaonPionIDnoSVD);
1068 variables->add(KaonProtonIDALL);
1069 variables->add(KaonProtonIDSVDonly);
1070 variables->add(KaonProtonIDnoSVD);
1071 variables->add(KaonElectronIDALL);
1072 variables->add(KaonElectronIDSVDonly);
1073 variables->add(KaonElectronIDnoSVD);
1074
1075 variables->add(KaonElectronLLSVDonly);
1076 variables->add(KaonPionLLSVDonly);
1077 variables->add(KaonKaonLLSVDonly);
1078 variables->add(KaonProtonLLSVDonly);
1079
1080 variables->add(KaonBinaryKaonPionIDALL);
1081 variables->add(KaonBinaryKaonPionIDSVDonly);
1082 variables->add(KaonBinaryKaonPionIDnoSVD);
1083 variables->add(KaonBinaryPionKaonIDALL);
1084 variables->add(KaonBinaryPionKaonIDSVDonly);
1085 variables->add(KaonBinaryPionKaonIDnoSVD);
1086 variables->add(KaonBinaryProtonKaonIDALL);
1087 variables->add(KaonBinaryProtonKaonIDSVDonly);
1088 variables->add(KaonBinaryProtonKaonIDnoSVD);
1089 variables->add(KaonBinaryElectronKaonIDALL);
1090 variables->add(KaonBinaryElectronKaonIDSVDonly);
1091 variables->add(KaonBinaryElectronKaonIDnoSVD);
1092
1093 variables->add(PionDPionIDALL);
1094 variables->add(PionDPionIDSVDonly);
1095 variables->add(PionDPionIDnoSVD);
1096 variables->add(PionDKaonIDALL);
1097 variables->add(PionDKaonIDSVDonly);
1098 variables->add(PionDKaonIDnoSVD);
1099 variables->add(PionDElectronIDALL);
1100 variables->add(PionDElectronIDSVDonly);
1101 variables->add(PionDElectronIDnoSVD);
1102 variables->add(PionDProtonIDALL);
1103 variables->add(PionDProtonIDSVDonly);
1104 variables->add(PionDProtonIDnoSVD);
1105
1106 variables->add(PionDElectronLLSVDonly);
1107 variables->add(PionDPionLLSVDonly);
1108 variables->add(PionDKaonLLSVDonly);
1109 variables->add(PionDProtonLLSVDonly);
1110
1111 variables->add(PionDBinaryPionKaonIDALL);
1112 variables->add(PionDBinaryPionKaonIDSVDonly);
1113 variables->add(PionDBinaryPionKaonIDnoSVD);
1114 variables->add(PionDBinaryKaonPionIDALL);
1115 variables->add(PionDBinaryKaonPionIDSVDonly);
1116 variables->add(PionDBinaryKaonPionIDnoSVD);
1117 variables->add(PionDBinaryProtonPionIDALL);
1118 variables->add(PionDBinaryProtonPionIDSVDonly);
1119 variables->add(PionDBinaryProtonPionIDnoSVD);
1120 variables->add(PionDBinaryElectronPionIDALL);
1121 variables->add(PionDBinaryElectronPionIDSVDonly);
1122 variables->add(PionDBinaryElectronPionIDnoSVD);
1123
1124 variables->add(SlowPionPionIDALL);
1125 variables->add(SlowPionPionIDSVDonly);
1126 variables->add(SlowPionPionIDnoSVD);
1127 variables->add(SlowPionKaonIDALL);
1128 variables->add(SlowPionKaonIDSVDonly);
1129 variables->add(SlowPionKaonIDnoSVD);
1130 variables->add(SlowPionElectronIDALL);
1131 variables->add(SlowPionElectronIDSVDonly);
1132 variables->add(SlowPionElectronIDnoSVD);
1133 variables->add(SlowPionProtonIDALL);
1134 variables->add(SlowPionProtonIDSVDonly);
1135 variables->add(SlowPionProtonIDnoSVD);
1136
1137 variables->add(SlowPionElectronLLSVDonly);
1138 variables->add(SlowPionPionLLSVDonly);
1139 variables->add(SlowPionKaonLLSVDonly);
1140 variables->add(SlowPionProtonLLSVDonly);
1141
1142 variables->add(SlowPionBinaryPionKaonIDALL);
1143 variables->add(SlowPionBinaryPionKaonIDSVDonly);
1144 variables->add(SlowPionBinaryPionKaonIDnoSVD);
1145 variables->add(SlowPionBinaryKaonPionIDALL);
1146 variables->add(SlowPionBinaryKaonPionIDSVDonly);
1147 variables->add(SlowPionBinaryKaonPionIDnoSVD);
1148 variables->add(SlowPionBinaryProtonPionIDALL);
1149 variables->add(SlowPionBinaryProtonPionIDSVDonly);
1150 variables->add(SlowPionBinaryProtonPionIDnoSVD);
1151 variables->add(SlowPionBinaryElectronPionIDALL);
1152 variables->add(SlowPionBinaryElectronPionIDSVDonly);
1153 variables->add(SlowPionBinaryElectronPionIDnoSVD);
1154
1155 RooDataSet* DstarDataset = new RooDataSet("DstarDataset", "DstarDataset", preselTree.get(), *variables);
1156
1157 if (DstarDataset->sumEntries() == 0) {
1158 B2FATAL("The Dstar dataset is empty, stopping here");
1159 }
1160
1161 RooPlot* DstarFitFrame = DstarDataset->plotOn(deltaM.frame());
1162
1163 RooRealVar GaussMean("GaussMean", "GaussMean", 0.145, 0.140, 0.150);
1164 RooRealVar GaussSigma1("GaussSigma1", "GaussSigma1", 0.01, 1.e-4, 1.0);
1165 RooGaussian DstarGauss1("DstarGauss1", "DstarGauss1", deltaM, GaussMean, GaussSigma1);
1166 RooRealVar GaussSigma2("GaussSigma2", "GaussSigma2", 0.001, 1.e-4, 1.0);
1167 RooGaussian DstarGauss2("DstarGauss2", "DstarGauss2", deltaM, GaussMean, GaussSigma2);
1168 RooRealVar fracGaussYield("fracGaussYield", "Fraction of two Gaussians", 0.75, 0.0, 1.0);
1169 RooAddPdf DstarSignalPDF("DstarSignalPDF", "DstarGauss1+DstarGauss2", RooArgList(DstarGauss1, DstarGauss2), fracGaussYield);
1170
1171 RooRealVar dm0Bkg("dm0Bkg", "dm0", 0.13957018, 0.130, 0.140);
1172 RooRealVar aBkg("aBkg", "a", -0.0784, -0.08, 3.0);
1173 RooRealVar bBkg("bBkg", "b", -0.444713, -0.5, 0.4);
1174 RooRealVar cBkg("cBkg", "c", 0.3);
1175 RooDstD0BG DstarBkgPDF("DstarBkgPDF", "DstarBkgPDF", deltaM, dm0Bkg, cBkg, aBkg, bBkg);
1176 RooRealVar nSignalDstar("nSignalDstar", "signal yield", 0.5 * preselTree->GetEntries(), 0, preselTree->GetEntries());
1177 RooRealVar nBkgDstar("nBkgDstar", "background yield", 0.5 * preselTree->GetEntries(), 0, preselTree->GetEntries());
1178 RooAddPdf totalPDFDstar("totalPDFDstar", "totalPDFDstar pdf", RooArgList(DstarSignalPDF, DstarBkgPDF),
1179 RooArgList(nSignalDstar, nBkgDstar));
1180
1181 B2INFO("Dstar: Start fitting...");
1182 RooFitResult* DstarFitResult = totalPDFDstar.fitTo(*DstarDataset, Save(kTRUE), PrintLevel(-1));
1183
1184 int status = DstarFitResult->status();
1185 int covqual = DstarFitResult->covQual();
1186 double diff = nSignalDstar.getValV() + nBkgDstar.getValV() - DstarDataset->sumEntries();
1187
1188 B2INFO("Dstar: Fit status: " << status << "; covariance quality: " << covqual);
1189 // if the fit is not healthy, try again once before giving up, with a slightly different setup:
1190 if ((status > 0) || (TMath::Abs(diff) > 1.) || (nSignalDstar.getError() < sqrt(nSignalDstar.getValV()))
1191 || (nSignalDstar.getError() > (nSignalDstar.getValV()))) {
1192
1193 DstarFitResult = totalPDFDstar.fitTo(*DstarDataset, Save(), Strategy(2), Offset(1));
1194 status = DstarFitResult->status();
1195 covqual = DstarFitResult->covQual();
1196 diff = nSignalDstar.getValV() + nBkgDstar.getValV() - DstarDataset->sumEntries();
1197 }
1198
1199 if ((status > 0) || (TMath::Abs(diff) > 1.) || (nSignalDstar.getError() < sqrt(nSignalDstar.getValV()))
1200 || (nSignalDstar.getError() > (nSignalDstar.getValV()))) {
1201 B2WARNING("Dstar: Fit problem: fit status " << status << "; sum of component yields minus the dataset yield is " << diff <<
1202 "; signal yield is " << nSignalDstar.getValV() << ", while its uncertainty is " << nSignalDstar.getError());
1203 }
1204 if (covqual < 2) {
1205 B2INFO("Dstar: Fit warning: covariance quality " << covqual);
1206 }
1207
1208 totalPDFDstar.plotOn(DstarFitFrame, LineColor(TColor::GetColor("#4575b4")));
1209
1210 double chisquare = DstarFitFrame->chiSquare();
1211 B2INFO("Dstar: Fit chi2 = " << chisquare);
1212 totalPDFDstar.paramOn(DstarFitFrame, Layout(0.63, 0.96, 0.93), Format("NEU", AutoPrecision(2)));
1213 DstarFitFrame->getAttText()->SetTextSize(0.03);
1214
1215 totalPDFDstar.plotOn(DstarFitFrame, Components("DstarSignalPDF"), LineColor(TColor::GetColor("#d73027")));
1216 totalPDFDstar.plotOn(DstarFitFrame, Components("DstarBkgPDF"), LineColor(TColor::GetColor("#fc8d59")));
1217 totalPDFDstar.plotOn(DstarFitFrame, LineColor(TColor::GetColor("#4575b4")));
1218
1219 DstarFitFrame->GetXaxis()->SetTitle("#Deltam [GeV/c^{2}]");
1220 TCanvas* canvDstar = new TCanvas("canvDstar", "canvDstar");
1221 canvDstar->cd();
1222
1223 DstarFitFrame->Draw();
1224
1225 if (m_isMakePlots) {
1226 canvDstar->Print("SVDdEdxValidationFitDstar.pdf");
1227 TFile DstarFitPlotFile("SVDdEdxValidationDstarFitPlotFile.root", "RECREATE");
1228 canvDstar->Write();
1229 DstarFitPlotFile.Close();
1230 }
1231
1233
1234 RooStats::SPlot* sPlotDatasetDstar = new RooStats::SPlot("sData", "An SPlot", *DstarDataset, &totalPDFDstar,
1235 RooArgList(nSignalDstar, nBkgDstar));
1236
1237 for (int iEvt = 0; iEvt < 5; iEvt++) {
1238 if (TMath::Abs(sPlotDatasetDstar->GetSWeight(iEvt, "nSignalDstar") + sPlotDatasetDstar->GetSWeight(iEvt, "nBkgDstar") - 1) > 5.e-3)
1239 B2FATAL("Dstar: sPlot error: sum of weights not equal to 1");
1240 }
1241
1242 RooDataSet* DstarDatasetSWeighted = new RooDataSet(DstarDataset->GetName(), DstarDataset->GetTitle(), DstarDataset,
1243 *DstarDataset->get());
1244
1245 TTree* treeDstar_sw = DstarDatasetSWeighted->GetClonedTree();
1246 treeDstar_sw->SetName("treeDstar_sw");
1247
1248 B2INFO("Dstar: sPlot done. ");
1249
1250 return treeDstar_sw;
1251}
Belle2::CalibrationAlgorithm
Base class for calibration algorithms.
Definition: CalibrationAlgorithm.h:37
Belle2::CalibrationAlgorithm::setDescription
void setDescription(const std::string &description)
Set algorithm description (in constructor)
Definition: CalibrationAlgorithm.h:321
Belle2::CalibrationAlgorithm::EResult
EResult
The result of calibration.
Definition: CalibrationAlgorithm.h:40
Belle2::CalibrationAlgorithm::c_OK
@ c_OK
Finished successfully =0 in Python.
Definition: CalibrationAlgorithm.h:41
Belle2::CalibrationAlgorithm::c_NotEnoughData
@ c_NotEnoughData
Needs more data =2 in Python.
Definition: CalibrationAlgorithm.h:43
Belle2::SVDdEdxValidationAlgorithm::m_MomLowROC
double m_MomLowROC
lower edge of the momentum interval considered for the ROC curve
Definition: SVDdEdxValidationAlgorithm.h:104
Belle2::SVDdEdxValidationAlgorithm::PlotROCCurve
void PlotROCCurve(TTree *SignalTree, TString SignalWeightName, TString SignalVarName, TString SignalVarNameFull, TTree *FakeTree, TString FakeWeightName, TString FakeVarName, TString FakeVarNameFull, TString PIDVarName)
a generic function to produce ROC curves
Definition: SVDdEdxValidationAlgorithm.cc:506
Belle2::SVDdEdxValidationAlgorithm::m_fullValidation
bool m_fullValidation
decide between full or basic validation mode.
Definition: SVDdEdxValidationAlgorithm.h:84
Belle2::SVDdEdxValidationAlgorithm::m_NumEffBins
unsigned int m_NumEffBins
number of bins for the efficiency/fake rate plot
Definition: SVDdEdxValidationAlgorithm.h:106
Belle2::SVDdEdxValidationAlgorithm::m_isMakePlots
bool m_isMakePlots
produce plots for monitoring of the fit quality
Definition: SVDdEdxValidationAlgorithm.h:85
Belle2::SVDdEdxValidationAlgorithm::m_MinEvtsPerTree
int m_MinEvtsPerTree
number of events in TTree below which we don't try to fit
Definition: SVDdEdxValidationAlgorithm.h:102
Belle2::SVDdEdxValidationAlgorithm::LambdaMassFit
TTree * LambdaMassFit(std::shared_ptr< TTree > preselTree)
produce histograms for protons
Definition: SVDdEdxValidationAlgorithm.cc:722
Belle2::SVDdEdxValidationAlgorithm::m_MomHighROC
double m_MomHighROC
upper edge of the momentum interval considered for the ROC curve
Definition: SVDdEdxValidationAlgorithm.h:105
Belle2::SVDdEdxValidationAlgorithm::m_MomHighEff
double m_MomHighEff
upper edge of the momentum interval for the efficiency/fake rate plot
Definition: SVDdEdxValidationAlgorithm.h:107
Belle2::SVDdEdxValidationAlgorithm::m_NumROCpoints
unsigned int m_NumROCpoints
number of points for the ROC curve plotting
Definition: SVDdEdxValidationAlgorithm.h:103
Belle2::SVDdEdxValidationAlgorithm::DstarMassFit
TTree * DstarMassFit(std::shared_ptr< TTree > preselTree)
produce histograms for K/pi(/mu)
Definition: SVDdEdxValidationAlgorithm.cc:922
Belle2::SVDdEdxValidationAlgorithm::calibrate
virtual EResult calibrate() override
run algorithm on data
Definition: SVDdEdxValidationAlgorithm.cc:54
Belle2::SVDdEdxValidationAlgorithm::PlotEfficiencyPlots
void PlotEfficiencyPlots(const TString &PIDDetectorsName, TTree *SignalTree, TString SignalWeightName, TString SignalVarName, TString SignalVarNameFull, TTree *FakeTree, TString FakeWeightName, TString FakeVarName, TString FakeVarNameFull, TString PIDVarName, TString PIDCut, unsigned int nbins, double MomLow, double MomHigh)
a generic function to produce efficiency plots
Definition: SVDdEdxValidationAlgorithm.cc:144
Belle2::sqrt
double sqrt(double a)
sqrt for double
Definition: beamHelpers.h:28
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:17