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
29#include <RooDataSet.h>
30#include <RooRealVar.h>
31#include <RooAddPdf.h>
32#include <RooGaussian.h>
33#include <RooChebychev.h>
34#include <RooBifurGauss.h>
35#include <RooDstD0BG.h>
36#include <RooAbsDataStore.h>
37#include <RooTreeDataStore.h>
38#include <RooMsgService.h>
39#include <RooStats/SPlot.h>
40
41using namespace RooFit;
42using namespace Belle2;
43
44SVDdEdxValidationAlgorithm::SVDdEdxValidationAlgorithm() : CalibrationAlgorithm("SVDdEdxValidationCollector"),
45 m_isMakePlots(true)
46{
47 setDescription("SVD dE/dx validation algorithm");
48}
49
50/* Main calibration method */
51CalibrationAlgorithm::EResult SVDdEdxValidationAlgorithm::calibrate()
52{
53 gROOT->SetBatch(true);
54
55 // const auto exprun = getRunList()[0];
56 // B2INFO("ExpRun used for calibration: " << exprun.first << " " << exprun.second);
57
58 // Get data objects
59 auto TTreeLambda = getObjectPtr<TTree>("Lambda");
60 auto TTreeDstar = getObjectPtr<TTree>("Dstar");
61 auto TTreeGamma = getObjectPtr<TTree>("Gamma");
62
63 if (TTreeLambda->GetEntries() < m_MinEvtsPerTree) {
64 B2WARNING("Not enough data for calibration.");
65 return c_NotEnoughData;
66 }
67
68 // call the calibration functions
69 TTree* TTreeLambdaSW = LambdaMassFit(TTreeLambda);
70 TTree* TTreeDstarSW = DstarMassFit(TTreeDstar);
71 TTree* TTreeGammaWrap = TTreeGamma.get();
72
73 std::vector<TString> PIDDetectors;
74 PIDDetectors.push_back("ALL");
75 PIDDetectors.push_back("SVDonly");
76 PIDDetectors.push_back("noSVD");
77
78 std::map<TTree*, TString> SWeightNameMap = {
79 {TTreeGammaWrap, "1"},
80 {TTreeDstarSW, "nSignalDstar_sw"},
81 {TTreeLambdaSW, "nSignalLambda_sw"}
82 };
83
84 for (const TString& PIDDetectorsName : PIDDetectors) {
85 PlotEfficiencyPlots(PIDDetectorsName, TTreeGammaWrap, SWeightNameMap[TTreeGammaWrap], "FirstElectron", "electron", TTreeDstarSW,
86 SWeightNameMap[TTreeDstarSW], "PionD", "pion",
87 "BinaryElectronPionID",
88 "0.5", m_NumEffBins, 0., m_MomHighEff);
89 PlotEfficiencyPlots(PIDDetectorsName, TTreeGammaWrap, SWeightNameMap[TTreeGammaWrap], "FirstElectron", "electron", TTreeDstarSW,
90 SWeightNameMap[TTreeDstarSW], "Kaon", "kaon",
91 "BinaryElectronKaonID", "0.5",
92 m_NumEffBins, 0., m_MomHighEff);
93 PlotEfficiencyPlots(PIDDetectorsName, TTreeLambdaSW, SWeightNameMap[TTreeLambdaSW], "Proton", "proton", TTreeDstarSW,
94 SWeightNameMap[TTreeDstarSW], "PionD", "pion",
95 "BinaryProtonPionID", "0.5",
96 m_NumEffBins, 0.25, m_MomHighEff);
97 PlotEfficiencyPlots(PIDDetectorsName, TTreeLambdaSW, SWeightNameMap[TTreeLambdaSW], "Proton", "proton", TTreeDstarSW,
98 SWeightNameMap[TTreeDstarSW], "Kaon", "kaon",
99 "BinaryProtonKaonID", "0.5",
100 m_NumEffBins, 0.25, m_MomHighEff);
101 PlotEfficiencyPlots(PIDDetectorsName, TTreeDstarSW, SWeightNameMap[TTreeDstarSW], "PionD", "pion", TTreeDstarSW,
102 SWeightNameMap[TTreeDstarSW],
103 "Kaon", "kaon",
104 "BinaryPionKaonID", "0.5", m_NumEffBins,
105 0., m_MomHighEff);
106 PlotEfficiencyPlots(PIDDetectorsName, TTreeDstarSW, SWeightNameMap[TTreeDstarSW], "Kaon", "kaon", TTreeDstarSW,
107 SWeightNameMap[TTreeDstarSW],
108 "PionD", "pion",
109 "BinaryKaonPionID", "0.5", m_NumEffBins,
110 0., m_MomHighEff);
111 }
112
113 PlotROCCurve(TTreeGammaWrap, SWeightNameMap[TTreeGammaWrap], "FirstElectron", "electron", TTreeDstarSW,
114 SWeightNameMap[TTreeDstarSW], "PionD",
115 "pion", "BinaryElectronPionID");
116 PlotROCCurve(TTreeGammaWrap, SWeightNameMap[TTreeGammaWrap], "FirstElectron", "electron", TTreeDstarSW,
117 SWeightNameMap[TTreeDstarSW], "Kaon",
118 "kaon", "BinaryElectronKaonID");
119 PlotROCCurve(TTreeLambdaSW, SWeightNameMap[TTreeLambdaSW], "Proton", "proton", TTreeDstarSW, SWeightNameMap[TTreeDstarSW], "PionD",
120 "pion",
121 "BinaryProtonPionID");
122 PlotROCCurve(TTreeLambdaSW, SWeightNameMap[TTreeLambdaSW], "Proton", "proton", TTreeDstarSW, SWeightNameMap[TTreeDstarSW], "Kaon",
123 "kaon",
124 "BinaryProtonKaonID");
125 PlotROCCurve(TTreeDstarSW, SWeightNameMap[TTreeDstarSW], "PionD", "pion", TTreeDstarSW, SWeightNameMap[TTreeDstarSW], "Kaon",
126 "kaon",
127 "BinaryPionKaonID");
128 PlotROCCurve(TTreeDstarSW, SWeightNameMap[TTreeDstarSW], "Kaon", "kaon", TTreeDstarSW, SWeightNameMap[TTreeDstarSW], "PionD",
129 "pion",
130 "BinaryKaonPionID");
131
132 B2INFO("SVD dE/dx validation done!");
133
134 return c_OK;
135}
136
137// generic efficiency and fake rate
138void SVDdEdxValidationAlgorithm::PlotEfficiencyPlots(const TString& PIDDetectorsName, TTree* SignalTree, TString SignalWeightName,
139 TString SignalVarName, TString SignalVarNameFull, TTree* FakeTree, TString FakeWeightName, TString FakeVarName,
140 TString FakeVarNameFull, TString PIDVarName, TString PIDCut, unsigned int nbins, double MomLow, double MomHigh)
141{
142
143 if ((SignalTree == nullptr) || (FakeTree == nullptr)) {
144 B2FATAL("Invalid dataset, stopping here");
145 }
146
147 if ((SignalTree->GetEntries() == 0) || (FakeTree->GetEntries() == 0)) {
148 B2FATAL("The dataset is empty, stopping here");
149 }
150
151 if ((SignalTree->GetBranch(Form("%sMomentum", SignalVarName.Data())) == nullptr)
152 || (FakeTree->GetBranch(Form("%sMomentum", FakeVarName.Data())) == nullptr)) {
153 B2FATAL("Check the provided branch name, stopping here");
154 }
155
156 TString SignalFiducialCut = "(1>0)"; // placeholder for a possible sanity cut
157 TString FakesFiducialCut = "(1>0)";
158
159 // Produce the plots of the SVD PID distribution
160 if (PIDDetectorsName == "SVDonly") {
161 SignalTree->Draw(Form("%s%s%s>>hSignalPIDDistribution(100,0.,1.)", SignalVarName.Data(), PIDVarName.Data(),
162 PIDDetectorsName.Data()),
163 SignalWeightName + Form("* (%sMomentum>%f && %sMomentum<%f)", SignalVarName.Data(), MomLow, SignalVarName.Data(), MomHigh), "goff");
164 TH1F* hSignalPIDDistribution = (TH1F*)gDirectory->Get("hSignalPIDDistribution");
165 hSignalPIDDistribution->Scale(1. / hSignalPIDDistribution->Integral());
166 hSignalPIDDistribution->GetXaxis()->SetTitle(PIDVarName + PIDDetectorsName + " for " + SignalVarNameFull);
167 hSignalPIDDistribution->GetYaxis()->SetTitle("Candidates, normalised");
168 hSignalPIDDistribution->SetMaximum(1.35 * hSignalPIDDistribution->GetMaximum());
169
170 TCanvas* DistribCanvas = new TCanvas("DistribCanvas", "", 600, 600);
171 gPad->SetTopMargin(0.05);
172 gPad->SetRightMargin(0.05);
173 gPad->SetLeftMargin(0.13);
174 gPad->SetBottomMargin(0.12);
175
176 hSignalPIDDistribution->SetLineWidth(2);
177 hSignalPIDDistribution->SetLineColor(TColor::GetColor("#2166ac"));
178 hSignalPIDDistribution->Draw("hist ");
179
180 DistribCanvas->Print("SVDdEdxValidation_Distribution_" + SignalVarNameFull + PIDVarName + PIDDetectorsName +
181 "_MomRange_" +
182 std::to_string(
183 MomLow)
184 .substr(0, 3) +
185 "_" + std::to_string(MomHigh).substr(0, 3) + ".pdf");
186 TFile DistribFile("SVDdEdxValidation_Distribution_" + SignalVarNameFull + PIDVarName + PIDDetectorsName +
187 "_MomRange_" +
188 std::to_string(
189 MomLow)
190 .substr(0, 3) +
191 "_" + std::to_string(MomHigh).substr(0, 3) + ".root", "RECREATE");
192 hSignalPIDDistribution->SetLineColor(kBlack);
193 hSignalPIDDistribution->Write();
194 DistribFile.Close();
195 delete DistribCanvas;
196 }
197
198 // ---------- Momentum distributions (for efficiency determination) ----------
199
200 SignalTree->Draw(Form("%sMomentum>>hAllSignal(%i,%f,%f)", SignalVarName.Data(), nbins, MomLow, MomHigh),
201 SignalWeightName + " * (" + SignalFiducialCut + ")", "goff");
202 SignalTree->Draw(Form("%sMomentum>>hSelectedSignal(%i,%f,%f)", SignalVarName.Data(), nbins, MomLow, MomHigh),
203 SignalWeightName + " * (" + SignalVarName + PIDVarName + PIDDetectorsName + ">" + PIDCut + "&&" + SignalFiducialCut +
204 ")",
205 "goff");
206
207 FakeTree->Draw(Form("%sMomentum>>hAllFakes(%i,%f,%f)", FakeVarName.Data(), nbins, MomLow, MomHigh),
208 FakeWeightName + " * (" + FakesFiducialCut + ")", "goff");
209 FakeTree->Draw(Form("%sMomentum>>hSelectedFakes(%i,%f,%f)", FakeVarName.Data(), nbins, MomLow, MomHigh),
210 FakeWeightName + " * (" + FakeVarName + PIDVarName + PIDDetectorsName + ">" + PIDCut + "&&" + FakesFiducialCut + ")",
211 "goff");
212
213 TH1F* hAllSignal = (TH1F*)gDirectory->Get("hAllSignal");
214 TH1F* hSelectedSignal = (TH1F*)gDirectory->Get("hSelectedSignal");
215 TH1F* hAllFakes = (TH1F*)gDirectory->Get("hAllFakes");
216 TH1F* hSelectedFakes = (TH1F*)gDirectory->Get("hSelectedFakes");
217
218 // ---------- Add slow pions to the pion dataset ----------
219 if (strncmp(SignalVarName.Data(), "PionD", 5) == 0) {
220 SignalTree->Draw(Form("SlowPionMomentum>>hAllSignalSlow(%i,%f,%f)", nbins, MomLow, MomHigh),
221 SignalWeightName + " * (" + SignalFiducialCut + ")", "goff");
222 SignalTree->Draw(Form("SlowPionMomentum>>hSelectedSignalSlow(%i,%f,%f)", nbins, MomLow, MomHigh),
223 SignalWeightName + " * (SlowPion" + PIDVarName + PIDDetectorsName + ">" + PIDCut + "&&" + SignalFiducialCut + ")", "goff");
224 TH1F* hAllSignalSlow = (TH1F*)gDirectory->Get("hAllSignalSlow");
225 TH1F* hSelectedSignalSlow = (TH1F*)gDirectory->Get("hSelectedSignalSlow");
226 hAllSignal->Add(hAllSignalSlow);
227 hSelectedSignal->Add(hSelectedSignalSlow);
228 }
229
230 if (strncmp(FakeVarName.Data(), "PionD", 5) == 0) {
231 FakeTree->Draw(Form("SlowPionMomentum>>hAllFakesSlow(%i,%f,%f)", nbins, MomLow, MomHigh),
232 FakeWeightName + " * (" + FakesFiducialCut + ")",
233 "goff");
234 FakeTree->Draw(Form("SlowPionMomentum>>hSelectedFakesSlow(%i,%f,%f)", nbins, MomLow, MomHigh),
235 FakeWeightName + " * (SlowPion" + PIDVarName + PIDDetectorsName + ">" + PIDCut + "&&" + FakesFiducialCut + ")", "goff");
236 TH1F* hAllFakesSlow = (TH1F*)gDirectory->Get("hAllFakesSlow");
237 TH1F* hSelectedFakesSlow = (TH1F*)gDirectory->Get("hSelectedFakesSlow");
238 hAllFakes->Add(hAllFakesSlow);
239 hSelectedFakes->Add(hSelectedFakesSlow);
240 }
241
242 TH1F* EffHistoSig = (TH1F*)hAllSignal->Clone("EffHistoSig"); // signal efficiency
243 TH1F* EffHistoFake = (TH1F*)hAllFakes->Clone("EffHistoFake"); // fakes efficiency
244
245 EffHistoSig->Divide(hSelectedSignal, hAllSignal, 1, 1, "B");
246 EffHistoFake->Divide(hSelectedFakes, hAllFakes, 1, 1, "B");
247
248 // PID plots
249 TH1F* hBase = new TH1F("hBase", "", 100, 0.0, MomHigh);
250 hBase->SetTitle(";Momentum [GeV];Efficiency");
251 hBase->SetMaximum(1.20);
252 hBase->SetMinimum(0.0);
253
254 TLegend* tleg1 = new TLegend(0.63, 0.82, 0.93, 0.94);
255 tleg1->AddEntry(EffHistoSig, SignalVarNameFull + " efficiency", "pl");
256 tleg1->AddEntry(EffHistoFake, FakeVarNameFull + " fake rate", "pl");
257
258 TCanvas* ResultCanvas = new TCanvas("ResultCanvas", "", 600, 600);
259 gPad->SetTopMargin(0.05);
260 gPad->SetRightMargin(0.05);
261 gPad->SetLeftMargin(0.13);
262 gPad->SetBottomMargin(0.12);
263
264 ResultCanvas->SetGrid();
265 hBase->Draw();
266 EffHistoSig->SetMarkerSize(1.5);
267 EffHistoSig->SetMarkerStyle(22);
268 EffHistoSig->SetMarkerColor(TColor::GetColor("#2166ac"));
269 EffHistoSig->SetLineColor(TColor::GetColor("#2166ac"));
270 EffHistoSig->Draw("P,same");
271
272 EffHistoFake->SetMarkerSize(1.5);
273 EffHistoFake->SetMarkerStyle(23);
274 EffHistoFake->SetMarkerColor(TColor::GetColor("#ef8a62"));
275 EffHistoFake->SetLineColor(TColor::GetColor("#ef8a62"));
276 EffHistoFake->Draw("P,same");
277
278 tleg1->Draw("same");
279
280 hBase->SetStats(0);
281 hBase->GetXaxis()->SetTitleSize(0.04);
282 hBase->GetYaxis()->SetTitleSize(0.04);
283 hBase->GetXaxis()->SetTitleOffset(1.0);
284 hBase->GetYaxis()->SetTitleOffset(1.3);
285 hBase->GetYaxis()->SetLabelSize(0.04);
286 hBase->GetXaxis()->SetLabelSize(0.04);
287
288 // std::setprecision(2);
289 ResultCanvas->Print("SVDdEdxValidation_Efficiency_" + SignalVarNameFull + "_vs_" + FakeVarNameFull + PIDVarName + "_" +
290 PIDDetectorsName +
291 "_Cut" +
292 PIDCut + "_MomRange_" + std::to_string(MomLow).substr(0, 3) + "_" + std::to_string(MomHigh).substr(0, 3) + ".pdf");
293 TFile ResultFile("SVDdEdxValidation_Efficiency_" + SignalVarNameFull + "_vs_" + FakeVarNameFull + PIDVarName + "_" +
294 PIDDetectorsName +
295 "_Cut" +
296 PIDCut + "_MomRange_" + std::to_string(MomLow).substr(0, 3) + "_" + std::to_string(MomHigh).substr(0, 3) + ".root", "RECREATE");
297 EffHistoSig->SetLineColor(kBlack);
298 EffHistoSig->SetMarkerColor(kBlack);
299 EffHistoFake->SetLineColor(kBlack);
300 EffHistoFake->SetMarkerColor(kBlack);
301 EffHistoSig->Write();
302 EffHistoFake->Write();
303 ResultFile.Close();
304 delete ResultCanvas;
305 delete hBase;
306}
307
308void SVDdEdxValidationAlgorithm::PlotROCCurve(TTree* SignalTree, TString SignalWeightName, TString SignalVarName,
309 TString SignalVarNameFull, TTree* FakeTree, TString FakeWeightName, TString FakeVarName, TString FakeVarNameFull,
310 TString PIDVarName)
311{
312
313 if ((SignalTree == nullptr) || (FakeTree == nullptr)) {
314 B2FATAL("Invalid dataset, stopping here");
315 }
316
317 if ((SignalTree->GetEntries() == 0) || (FakeTree->GetEntries() == 0)) {
318 B2FATAL("The dataset is empty, stopping here");
319 }
320
321 if ((SignalTree->GetBranch(Form("%sMomentum", SignalVarName.Data())) == nullptr)
322 || (FakeTree->GetBranch(Form("%sMomentum", FakeVarName.Data())) == nullptr)) {
323 B2FATAL("Check the provided branch name, stopping here");
324 }
325
326 std::vector<TString> PIDDetectors;
327 PIDDetectors.clear();
328 PIDDetectors.push_back("ALL");
329 PIDDetectors.push_back("noSVD");
330
331 std::vector<double> SignalEfficiencyALL, FakeEfficiencyALL;
332 SignalEfficiencyALL.reserve(m_NumROCpoints);
333 FakeEfficiencyALL.reserve(m_NumROCpoints);
334 std::vector<double> SignalEfficiencynoSVD, FakeEfficiencynoSVD;
335 SignalEfficiencynoSVD.reserve(m_NumROCpoints);
336 FakeEfficiencynoSVD.reserve(m_NumROCpoints);
337
338 TString SignalFiducialCut = SignalVarName + PIDVarName + "noSVD>=0"; // sanity cuts to reject events with NaN
339 TString FakesFiducialCut = FakeVarName + PIDVarName + "noSVD>=0";
340
341 // calculate efficiencies
342 for (unsigned int i = 0; i < PIDDetectors.size(); i++) {
343 for (unsigned int j = 0; j < m_NumROCpoints; ++j) {
344 delete gROOT->FindObject("PIDCut");
345 delete gROOT->FindObject("hAllSignal");
346 delete gROOT->FindObject("hSelectedSignal");
347
348 // scan cut values from 0 to 1, with a denser scan closer to 0 or 1, to get a nicer ROC curve
349 double x = 1. / m_NumROCpoints * j;
350 TString PIDCut = TString::Format("%f", 1. / (1 + TMath::Power(x / (1 - x), -3)));
351
352 // TString PIDCut = TString::Format("%f", 0. + 1. / m_NumROCpoints * j);
353
354 SignalTree->Draw(Form("%sMomentum>>hAllSignal(1,%f,%f)", SignalVarName.Data(), m_MomLowROC, m_MomHighROC),
355 SignalWeightName + " * (" + SignalFiducialCut + ")", "goff");
356 SignalTree->Draw(Form("%sMomentum>>hSelectedSignal(1,%f,%f)", SignalVarName.Data(), m_MomLowROC, m_MomHighROC),
357 SignalWeightName + " * (" + SignalVarName + PIDVarName + PIDDetectors[i] + ">" + PIDCut + "&&" + SignalFiducialCut + ")",
358 "goff");
359
360 TH1F* hAllSignal = (TH1F*)gDirectory->Get("hAllSignal");
361 TH1F* hSelectedSignal = (TH1F*)gDirectory->Get("hSelectedSignal");
362
363 if (strncmp(SignalVarName.Data(), "PionD", 5) == 0) {
364 SignalTree->Draw(Form("SlowPionMomentum>>hAllSignalSlow(1,%f,%f)", m_MomLowROC, m_MomHighROC),
365 SignalWeightName + " * (" + SignalFiducialCut + "&& SlowPion" + PIDVarName + "noSVD>=0" + ")", "goff");
366 SignalTree->Draw(Form("SlowPionMomentum>>hSelectedSignalSlow(1,%f,%f)", m_MomLowROC, m_MomHighROC),
367 SignalWeightName + " * (SlowPion" + PIDVarName + PIDDetectors[i] + ">" + PIDCut + "&&" + SignalFiducialCut + "&& SlowPion" +
368 PIDVarName +
369 "noSVD>=0" + ")", "goff");
370 TH1F* hAllSignalSlow = (TH1F*)gDirectory->Get("hAllSignalSlow");
371 TH1F* hSelectedSignalSlow = (TH1F*)gDirectory->Get("hSelectedSignalSlow");
372 hAllSignal->Add(hAllSignalSlow);
373 hSelectedSignal->Add(hSelectedSignalSlow);
374 }
375
376 if (PIDDetectors[i] == "ALL") {
377 SignalEfficiencyALL.push_back(hSelectedSignal->Integral() / hAllSignal->Integral());
378 }
379
380 if (PIDDetectors[i] == "noSVD") {
381 SignalEfficiencynoSVD.push_back(hSelectedSignal->Integral() / hAllSignal->Integral());
382 }
383 }
384 }
385
386 // calculate fake rates
387
388 for (unsigned int i = 0; i < PIDDetectors.size(); i++) {
389 for (unsigned int j = 0; j < m_NumROCpoints; ++j) {
390 delete gROOT->FindObject("PIDCut");
391 delete gROOT->FindObject("hAllFakes");
392 delete gROOT->FindObject("hSelectedFakes");
393
394 // scan cut values from 0 to 1, with a denser scan closer to 0 or 1, to get a nicer ROC curve
395 double x = 1. / m_NumROCpoints * j;
396 TString PIDCut = TString::Format("%f", 1. / (1 + TMath::Power(x / (1 - x), -3)));
397
398 FakeTree->Draw(Form("%sMomentum>>hAllFakes(1,%f,%f)", FakeVarName.Data(), m_MomLowROC, m_MomHighROC),
399 FakeWeightName + " * (" + FakesFiducialCut + ")", "goff");
400 FakeTree->Draw(Form("%sMomentum>>hSelectedFakes(1,%f,%f)", FakeVarName.Data(), m_MomLowROC, m_MomHighROC),
401 FakeWeightName + " * (" + FakeVarName + PIDVarName + PIDDetectors[i] + ">" + PIDCut + "&&" + FakesFiducialCut + ")", "goff");
402
403 TH1F* hSelectedFakes = (TH1F*)gDirectory->Get("hSelectedFakes");
404 TH1F* hAllFakes = (TH1F*)gDirectory->Get("hAllFakes");
405
406 if (strncmp(FakeVarName.Data(), "PionD", 5) == 0) {
407 FakeTree->Draw(Form("SlowPionMomentum>>hAllFakesSlow(1,%f,%f)", m_MomLowROC, m_MomHighROC),
408 FakeWeightName + " * (" + FakesFiducialCut + "&& SlowPion" + PIDVarName + "noSVD>=0" + ")", "goff");
409 FakeTree->Draw(Form("SlowPionMomentum>>hSelectedFakesSlow(1,%f,%f)", m_MomLowROC, m_MomHighROC),
410 FakeWeightName + " * (SlowPion" + PIDVarName + PIDDetectors[i] + ">" + PIDCut + "&&" + FakesFiducialCut + "&& SlowPion" + PIDVarName
411 +
412 "noSVD>=0" + ")", "goff");
413 TH1F* hAllFakesSlow = (TH1F*)gDirectory->Get("hAllFakesSlow");
414 TH1F* hSelectedFakesSlow = (TH1F*)gDirectory->Get("hSelectedFakesSlow");
415 hAllFakes->Add(hAllFakesSlow);
416 hSelectedFakes->Add(hSelectedFakesSlow);
417 }
418
419 if (PIDDetectors[i] == "ALL") {
420 FakeEfficiencyALL.push_back(hSelectedFakes->Integral() / hAllFakes->Integral());
421 }
422
423 if (PIDDetectors[i] == "noSVD") {
424 FakeEfficiencynoSVD.push_back(hSelectedFakes->Integral() / hAllFakes->Integral());
425 }
426 }
427 }
428
429 auto ResultCanvas = new TCanvas("ResultCanvas", "", 600, 400);
430 TMultiGraph* hmgraph = new TMultiGraph();
431
432 // efficiency and kaon fake rate
433 TGraph* hgraphALL = new TGraph(m_NumROCpoints, FakeEfficiencyALL.data(), SignalEfficiencyALL.data());
434 hgraphALL->SetMarkerColor(TColor::GetColor("#2166ac"));
435 hgraphALL->SetMarkerStyle(20);
436 hgraphALL->SetLineColor(TColor::GetColor("#2166ac"));
437 hgraphALL->SetLineWidth(3);
438 hgraphALL->SetDrawOption("AP*");
439 hgraphALL->SetTitle("with SVD");
440
441 TGraph* hgraphnoSVD = new TGraph(m_NumROCpoints, FakeEfficiencynoSVD.data(), SignalEfficiencynoSVD.data());
442 hgraphnoSVD->SetMarkerColor(TColor::GetColor("#ef8a62"));
443 hgraphnoSVD->SetLineColor(TColor::GetColor("#ef8a62"));
444 hgraphnoSVD->SetLineWidth(3);
445 hgraphnoSVD->SetMarkerStyle(22);
446 hgraphnoSVD->SetDrawOption("P*");
447 hgraphnoSVD->SetTitle("without SVD");
448
449 hmgraph->Add(hgraphALL);
450 hmgraph->Add(hgraphnoSVD);
451 hmgraph->Draw("A");
452 hmgraph->GetHistogram()->GetXaxis()->SetTitle(FakeVarNameFull + " fake rate");
453 hmgraph->GetHistogram()->GetYaxis()->SetTitle(SignalVarNameFull + " signal efficiency");
454
455 ResultCanvas->BuildLegend(0.6, 0.25, 0.9, 0.5);
456 ResultCanvas->SetGrid();
457
458 ResultCanvas->Print("SVDdEdxValidation_ROC_curve_" + SignalVarNameFull + "_vs_" + FakeVarNameFull + PIDVarName + "_MomRange" +
459 std::to_string(m_MomLowROC).substr(0, 3) + "_" + std::to_string(m_MomHighROC).substr(0, 3) + ".pdf");
460
461 TFile ResultFile("SVDdEdxValidation_ROC_curve_" + SignalVarNameFull + "_vs_" + FakeVarNameFull + PIDVarName + "_MomRange" +
462 std::to_string(m_MomLowROC).substr(0, 3) + "_" + std::to_string(m_MomHighROC).substr(0, 3) + ".root", "RECREATE");
463 hmgraph->Write();
464 ResultFile.Close();
465
466 delete ResultCanvas;
467}
468
469TTree* SVDdEdxValidationAlgorithm::LambdaMassFit(std::shared_ptr<TTree> preselTree)
470{
471 B2INFO("Configuring the Lambda fit...");
472 gROOT->SetBatch(true);
473 RooMsgService::instance().setGlobalKillBelow(RooFit::WARNING);
474
475 RooRealVar InvM("InvM", "m(p^{+}#pi^{-})", 1.1, 1.13, "GeV/c^{2}");
476
477 RooRealVar ProtonMomentum("ProtonMomentum", "momentum for p", -1.e8, 1.e8);
478 RooRealVar ProtonSVDdEdx("ProtonSVDdEdx", "", -1.e8, 1.e8);
479
480 RooRealVar exp("exp", "experiment number", 0, 1.e5);
481 RooRealVar run("run", "run number", 0, 1.e7);
482
483 RooRealVar ProtonProtonIDALL("ProtonProtonIDALL", "", -1.e8, 1.e8);
484 RooRealVar ProtonProtonIDSVDonly("ProtonProtonIDSVDonly", "", -1.e8, 1.e8);
485 RooRealVar ProtonProtonIDnoSVD("ProtonProtonIDnoSVD", "", -1.e8, 1.e8);
486
487 RooRealVar ProtonBinaryProtonPionIDALL("ProtonBinaryProtonPionIDALL", "", -1.e8, 1.e8);
488 RooRealVar ProtonBinaryProtonPionIDSVDonly("ProtonBinaryProtonPionIDSVDonly", "", -1.e8, 1.e8);
489 RooRealVar ProtonBinaryProtonPionIDnoSVD("ProtonBinaryProtonPionIDnoSVD", "", -1.e8, 1.e8);
490
491 RooRealVar ProtonBinaryProtonKaonIDALL("ProtonBinaryProtonKaonIDALL", "", -1.e8, 1.e8);
492 RooRealVar ProtonBinaryProtonKaonIDSVDonly("ProtonBinaryProtonKaonIDSVDonly", "", -1.e8, 1.e8);
493 RooRealVar ProtonBinaryProtonKaonIDnoSVD("ProtonBinaryProtonKaonIDnoSVD", "", -1.e8, 1.e8);
494
495 RooRealVar ProtonBinaryProtonElectronIDALL("ProtonBinaryProtonElectronIDALL", "", -1.e8, 1.e8);
496 RooRealVar ProtonBinaryProtonElectronIDSVDonly("ProtonBinaryProtonElectronIDSVDonly", "", -1.e8, 1.e8);
497 RooRealVar ProtonBinaryProtonElectronIDnoSVD("ProtonBinaryProtonElectronIDnoSVD", "", -1.e8, 1.e8);
498
499 RooRealVar ProtonBinaryPionProtonIDALL("ProtonBinaryPionProtonIDALL", "", -1.e8, 1.e8);
500 RooRealVar ProtonBinaryPionProtonIDSVDonly("ProtonBinaryPionProtonIDSVDonly", "", -1.e8, 1.e8);
501 RooRealVar ProtonBinaryPionProtonIDnoSVD("ProtonBinaryPionProtonIDnoSVD", "", -1.e8, 1.e8);
502
503 RooRealVar ProtonBinaryKaonProtonIDALL("ProtonBinaryKaonProtonIDALL", "", -1.e8, 1.e8);
504 RooRealVar ProtonBinaryKaonProtonIDSVDonly("ProtonBinaryKaonProtonIDSVDonly", "", -1.e8, 1.e8);
505 RooRealVar ProtonBinaryKaonProtonIDnoSVD("ProtonBinaryKaonProtonIDnoSVD", "", -1.e8, 1.e8);
506
507 RooRealVar ProtonBinaryElectronProtonIDALL("ProtonBinaryElectronProtonIDALL", "", -1.e8, 1.e8);
508 RooRealVar ProtonBinaryElectronProtonIDSVDonly("ProtonBinaryElectronProtonIDSVDonly", "", -1.e8, 1.e8);
509 RooRealVar ProtonBinaryElectronProtonIDnoSVD("ProtonBinaryElectronProtonIDnoSVD", "", -1.e8, 1.e8);
510
511 auto variables = new RooArgSet();
512
513 variables->add(InvM);
514
515 variables->add(ProtonMomentum);
516 variables->add(ProtonSVDdEdx);
517 variables->add(exp);
518 variables->add(run);
519
520 variables->add(ProtonProtonIDALL);
521 variables->add(ProtonProtonIDSVDonly);
522 variables->add(ProtonProtonIDnoSVD);
523 variables->add(ProtonBinaryProtonPionIDALL);
524 variables->add(ProtonBinaryProtonPionIDSVDonly);
525 variables->add(ProtonBinaryProtonPionIDnoSVD);
526 variables->add(ProtonBinaryProtonKaonIDALL);
527 variables->add(ProtonBinaryProtonKaonIDSVDonly);
528 variables->add(ProtonBinaryProtonKaonIDnoSVD);
529 variables->add(ProtonBinaryProtonElectronIDALL);
530 variables->add(ProtonBinaryProtonElectronIDSVDonly);
531 variables->add(ProtonBinaryProtonElectronIDnoSVD);
532 variables->add(ProtonBinaryPionProtonIDALL);
533 variables->add(ProtonBinaryPionProtonIDSVDonly);
534 variables->add(ProtonBinaryPionProtonIDnoSVD);
535 variables->add(ProtonBinaryKaonProtonIDALL);
536 variables->add(ProtonBinaryKaonProtonIDSVDonly);
537 variables->add(ProtonBinaryKaonProtonIDnoSVD);
538 variables->add(ProtonBinaryElectronProtonIDALL);
539 variables->add(ProtonBinaryElectronProtonIDSVDonly);
540 variables->add(ProtonBinaryElectronProtonIDnoSVD);
541
542 RooDataSet* LambdaDataset = new RooDataSet("LambdaDataset", "LambdaDataset", preselTree.get(), *variables);
543
544 if (LambdaDataset->sumEntries() == 0) {
545 B2FATAL("The Lambda dataset is empty, stopping here");
546 }
547
548 // the signal PDF; might be revisited at a later point
549
550 RooRealVar GaussMean("GaussMean", " GaussMean", 1.116, 1.111, 1.12);
551 RooRealVar GaussSigma("GaussSigma", "#sigma_{1}", 3.e-3, 3.e-5, 10.e-3);
552 RooGaussian LambdaGauss("LambdaGauss", "LambdaGauss", InvM, GaussMean, GaussSigma);
553
554 /* temporary RooRealVar sigmaBifurGaussL1 and sigmaBifurGaussR1 to replace
555 * RooRealVar resolutionParamL("resolutionParamL", "resolutionParamL", 0.4, 5.e-4, 1.0);
556 * RooRealVar resolutionParamR("resolutionParamR", "resolutionParamR", 0.4, 5.e-4, 1.0);
557 * RooFormulaVar sigmaBifurGaussL1("sigmaBifurGaussL1", "resolutionParamL*GaussSigma", RooArgSet(resolutionParamL, GaussSigma));
558 * RooFormulaVar sigmaBifurGaussR1("sigmaBifurGaussR1", "resolutionParamR*GaussSigma", RooArgSet(resolutionParamR, GaussSigma));
559 */
560 RooRealVar sigmaBifurGaussL1("sigmaBifurGaussL1", "sigma left", 0.4 * 3.e-3, 3.e-5, 10.e-3);
561 RooRealVar sigmaBifurGaussR1("sigmaBifurGaussR1", "sigma right", 0.4 * 3.e-3, 3.e-5, 10.e-3);
562 RooBifurGauss LambdaBifurGauss("LambdaBifurGauss", "LambdaBifurGauss", InvM, GaussMean, sigmaBifurGaussL1, sigmaBifurGaussR1);
563
564 /* temporary RooRealVar sigmaBifurGaussL2 to replace
565 * RooRealVar resolutionParam2("resolutionParam2", "resolutionParam2", 0.2, 5.e-4, 1.0);
566 * sigmaBifurGaussL2("sigmaBifurGaussL2", "resolutionParam2*GaussSigma", RooArgSet(resolutionParam2, GaussSigma));
567 */
568 RooRealVar sigmaBifurGaussL2("sigmaBifurGaussL2", "sigmaBifurGaussL2", 0.2 * 3.e-3, 3.e-5, 10.e-3);
569 RooGaussian LambdaBifurGauss2("LambdaBifurGauss2", "LambdaBifurGauss2", InvM, GaussMean, sigmaBifurGaussL2);
570
571 RooRealVar fracBifurGaussYield("fracBifurGaussYield", "fracBifurGaussYield", 0.3, 5.e-4, 1.0);
572 RooRealVar fracGaussYield("fracGaussYield", "fracGaussYield", 0.8, 5.e-4, 1.0);
573
574 RooAddPdf LambdaCombinedBifurGauss("LambdaCombinedBifurGauss", "LambdaBifurGauss + LambdaBifurGauss2 ", RooArgList(LambdaBifurGauss,
575 LambdaBifurGauss2), RooArgList(fracBifurGaussYield));
576
577 RooAddPdf LambdaSignalPDF("LambdaSignalPDF", "LambdaCombinedBifurGauss + LambdaGauss", RooArgList(LambdaCombinedBifurGauss,
578 LambdaGauss), RooArgList(fracGaussYield));
579
580 // Background PDF
581 RooRealVar BkgPolyCoef0("BkgPolyCoef0", "BkgPolyCoef0", 0.1, 0., 1.5);
582 RooRealVar BkgPolyCoef1("BkgPolyCoef1", "BkgPolyCoef1", -0.5, -1.5, -1.e-3);
583 RooChebychev BkgPolyPDF("BkgPolyPDF", "BkgPolyPDF", InvM, RooArgList(BkgPolyCoef0, BkgPolyCoef1));
584
585 RooRealVar nSignalLambda("nSignalLambda", "nSignalLambda", 0.6 * preselTree->GetEntries(), 0., 0.99 * preselTree->GetEntries());
586 RooRealVar nBkgLambda("nBkgLambda", "nBkgLambda", 0.4 * preselTree->GetEntries(), 0., 0.99 * preselTree->GetEntries());
587 RooAddPdf totalPDFLambda("totalPDFLambda", "totalPDFLambda pdf", RooArgList(LambdaSignalPDF, BkgPolyPDF),
588 RooArgList(nSignalLambda, nBkgLambda));
589
590 B2INFO("Lambda: Start fitting...");
591 RooFitResult* LambdaFitResult = totalPDFLambda.fitTo(*LambdaDataset, Save(kTRUE), PrintLevel(-1));
592
593 int status = LambdaFitResult->status();
594 int covqual = LambdaFitResult->covQual();
595 double diff = nSignalLambda.getValV() + nBkgLambda.getValV() - LambdaDataset->sumEntries();
596
597 B2INFO("Lambda: Fit status: " << status << "; covariance quality: " << covqual);
598 // if the fit is not healthy, try again once before giving up, with a slightly different setup:
599 if ((status > 0) || (TMath::Abs(diff) > 1.) || (nSignalLambda.getError() < sqrt(nSignalLambda.getValV()))
600 || (nSignalLambda.getError() > (nSignalLambda.getValV()))) {
601
602 LambdaFitResult = totalPDFLambda.fitTo(*LambdaDataset, Save(), Strategy(2), Offset(1));
603 status = LambdaFitResult->status();
604 covqual = LambdaFitResult->covQual();
605 diff = nSignalLambda.getValV() + nBkgLambda.getValV() - LambdaDataset->sumEntries();
606 }
607
608 if ((status > 0) || (TMath::Abs(diff) > 1.) || (nSignalLambda.getError() < sqrt(nSignalLambda.getValV()))
609 || (nSignalLambda.getError() > (nSignalLambda.getValV()))) {
610 B2WARNING("Lambda: Fit problem: fit status " << status << "; sum of component yields minus the dataset yield is " << diff <<
611 "; signal yield is " << nSignalLambda.getValV() << ", while its uncertainty is " << nSignalLambda.getError());
612 }
613 if (covqual < 2) {
614 B2INFO("Lambda: Fit warning: covariance quality " << covqual);
615 }
616
617 TCanvas* canvLambda = new TCanvas("canvLambda", "canvLambda");
618 RooPlot* LambdaFitFrame = LambdaDataset->plotOn(InvM.frame(130));
619 totalPDFLambda.plotOn(LambdaFitFrame, LineColor(TColor::GetColor("#4575b4")));
620
621 double chisquare = LambdaFitFrame->chiSquare();
622 B2INFO("Lambda: Fit chi2 = " << chisquare);
623 totalPDFLambda.paramOn(LambdaFitFrame, Layout(0.6, 0.96, 0.93), Format("NEU", AutoPrecision(2)));
624 LambdaFitFrame->getAttText()->SetTextSize(0.03);
625
626 totalPDFLambda.plotOn(LambdaFitFrame, Components("LambdaSignalPDF"), LineColor(TColor::GetColor("#d73027")));
627 totalPDFLambda.plotOn(LambdaFitFrame, Components("BkgPolyPDF"), LineColor(TColor::GetColor("#fc8d59")));
628 totalPDFLambda.plotOn(LambdaFitFrame, LineColor(TColor::GetColor("#4575b4")));
629
630 LambdaFitFrame->GetXaxis()->SetTitle("m(p#pi^{-}) (GeV/c^{2})");
631
632 LambdaFitFrame->Draw();
633
634 if (m_isMakePlots) {
635 canvLambda->Print("SVDdEdxValidationFitLambda.pdf");
636 TFile LambdaFitPlotFile("SVDdEdxValidationLambdaFitPlotFile.root", "RECREATE");
637 canvLambda->Write();
638 LambdaFitPlotFile.Close();
639 }
640 RooStats::SPlot* sPlotDatasetLambda = new RooStats::SPlot("sData", "An SPlot", *LambdaDataset, &totalPDFLambda,
641 RooArgList(nSignalLambda, nBkgLambda));
642
643 for (int iEvt = 0; iEvt < 5; iEvt++) {
644 if (TMath::Abs(sPlotDatasetLambda->GetSWeight(iEvt, "nSignalLambda") + sPlotDatasetLambda->GetSWeight(iEvt,
645 "nBkgLambda") - 1) > 5.e-3)
646 B2FATAL("Lambda: sPlot error: sum of weights not equal to 1");
647 }
648
649 RooDataSet* LambdaDatasetSWeighted = new RooDataSet(LambdaDataset->GetName(), LambdaDataset->GetTitle(), LambdaDataset,
650 *LambdaDataset->get());
651
652 RooDataSet::setDefaultStorageType(RooAbsData::Tree);
653 ((RooTreeDataStore*)(LambdaDatasetSWeighted->store())->tree())->SetName("treeLambda_sw");
654 TTree* treeLambda_sw = LambdaDatasetSWeighted->GetClonedTree();
655
656 B2INFO("Lambda: sPlot done. ");
657
658 return treeLambda_sw;
659}
660
661TTree* SVDdEdxValidationAlgorithm::DstarMassFit(std::shared_ptr<TTree> preselTree)
662{
663 B2INFO("Configuring the Dstar fit...");
664 gROOT->SetBatch(true);
665 RooMsgService::instance().setGlobalKillBelow(RooFit::WARNING);
666
667 RooRealVar deltaM("deltaM", "m(D*)-m(D^{0})", 0.139545, 0.151, "GeV/c^{2}");
668
669 RooRealVar KaonMomentum("KaonMomentum", "momentum for Kaon(GeV)", -1.e8, 1.e8);
670 RooRealVar KaonSVDdEdx("KaonSVDdEdx", "", -1.e8, 1.e8);
671 RooRealVar PionDMomentum("PionDMomentum", "momentum for pion(GeV)", -1.e8, 1.e8);
672 RooRealVar PionDSVDdEdx("PionDSVDdEdx", "", -1.e8, 1.e8);
673 RooRealVar SlowPionMomentum("SlowPionMomentum", "momentum for slow pion(GeV)", -1.e8, 1.e8);
674 RooRealVar SlowPionSVDdEdx("SlowPionSVDdEdx", "", -1.e8, 1.e8);
675
676 RooRealVar exp("exp", "experiment number", 0, 1.e5);
677 RooRealVar run("run", "run number", 0, 1.e8);
678
679 RooRealVar KaonKaonIDALL("KaonKaonIDALL", "", -1.e8, 1.e8);
680 RooRealVar KaonKaonIDSVDonly("KaonKaonIDSVDonly", "", -1.e8, 1.e8);
681 RooRealVar KaonKaonIDnoSVD("KaonKaonIDnoSVD", "", -1.e8, 1.e8);
682
683 RooRealVar KaonPionIDALL("KaonPionIDALL", "", -1.e8, 1.e8);
684 RooRealVar KaonPionIDSVDonly("KaonPionIDSVDonly", "", -1.e8, 1.e8);
685 RooRealVar KaonPionIDnoSVD("KaonPionIDnoSVD", "", -1.e8, 1.e8);
686
687 RooRealVar KaonProtonIDALL("KaonProtonIDALL", "", -1.e8, 1.e8);
688 RooRealVar KaonProtonIDSVDonly("KaonProtonIDSVDonly", "", -1.e8, 1.e8);
689 RooRealVar KaonProtonIDnoSVD("KaonProtonIDnoSVD", "", -1.e8, 1.e8);
690
691 RooRealVar KaonElectronIDALL("KaonElectronIDALL", "", -1.e8, 1.e8);
692 RooRealVar KaonElectronIDSVDonly("KaonElectronIDSVDonly", "", -1.e8, 1.e8);
693 RooRealVar KaonElectronIDnoSVD("KaonElectronIDnoSVD", "", -1.e8, 1.e8);
694
695 RooRealVar KaonBinaryKaonPionIDALL("KaonBinaryKaonPionIDALL", "", -1.e8, 1.e8);
696 RooRealVar KaonBinaryKaonPionIDSVDonly("KaonBinaryKaonPionIDSVDonly", "", -1.e8, 1.e8);
697 RooRealVar KaonBinaryKaonPionIDnoSVD("KaonBinaryKaonPionIDnoSVD", "", -1.e8, 1.e8);
698
699 RooRealVar KaonBinaryPionKaonIDALL("KaonBinaryPionKaonIDALL", "", -1.e8, 1.e8);
700 RooRealVar KaonBinaryPionKaonIDSVDonly("KaonBinaryPionKaonIDSVDonly", "", -1.e8, 1.e8);
701 RooRealVar KaonBinaryPionKaonIDnoSVD("KaonBinaryPionKaonIDnoSVD", "", -1.e8, 1.e8);
702
703 RooRealVar KaonBinaryProtonKaonIDALL("KaonBinaryProtonKaonIDALL", "", -1.e8, 1.e8);
704 RooRealVar KaonBinaryProtonKaonIDSVDonly("KaonBinaryProtonKaonIDSVDonly", "", -1.e8, 1.e8);
705 RooRealVar KaonBinaryProtonKaonIDnoSVD("KaonBinaryProtonKaonIDnoSVD", "", -1.e8, 1.e8);
706
707 RooRealVar KaonBinaryElectronKaonIDALL("KaonBinaryElectronKaonIDALL", "", -1.e8, 1.e8);
708 RooRealVar KaonBinaryElectronKaonIDSVDonly("KaonBinaryElectronKaonIDSVDonly", "", -1.e8, 1.e8);
709 RooRealVar KaonBinaryElectronKaonIDnoSVD("KaonBinaryElectronKaonIDnoSVD", "", -1.e8, 1.e8);
710
711 RooRealVar PionDKaonIDALL("PionDKaonIDALL", "", -1.e8, 1.e8);
712 RooRealVar PionDKaonIDSVDonly("PionDKaonIDSVDonly", "", -1.e8, 1.e8);
713 RooRealVar PionDKaonIDnoSVD("PionDKaonIDnoSVD", "", -1.e8, 1.e8);
714
715 RooRealVar PionDPionIDALL("PionDPionIDALL", "", -1.e8, 1.e8);
716 RooRealVar PionDPionIDSVDonly("PionDPionIDSVDonly", "", -1.e8, 1.e8);
717 RooRealVar PionDPionIDnoSVD("PionDPionIDnoSVD", "", -1.e8, 1.e8);
718
719 RooRealVar PionDElectronIDALL("PionDElectronIDALL", "", -1.e8, 1.e8);
720 RooRealVar PionDElectronIDSVDonly("PionDElectronIDSVDonly", "", -1.e8, 1.e8);
721 RooRealVar PionDElectronIDnoSVD("PionDElectronIDnoSVD", "", -1.e8, 1.e8);
722
723 RooRealVar PionDProtonIDALL("PionDProtonIDALL", "", -1.e8, 1.e8);
724 RooRealVar PionDProtonIDSVDonly("PionDProtonIDSVDonly", "", -1.e8, 1.e8);
725 RooRealVar PionDProtonIDnoSVD("PionDProtonIDnoSVD", "", -1.e8, 1.e8);
726
727 RooRealVar PionDBinaryPionKaonIDALL("PionDBinaryPionKaonIDALL", "", -1.e8, 1.e8);
728 RooRealVar PionDBinaryPionKaonIDSVDonly("PionDBinaryPionKaonIDSVDonly", "", -1.e8, 1.e8);
729 RooRealVar PionDBinaryPionKaonIDnoSVD("PionDBinaryPionKaonIDnoSVD", "", -1.e8, 1.e8);
730
731 RooRealVar PionDBinaryKaonPionIDALL("PionDBinaryKaonPionIDALL", "", -1.e8, 1.e8);
732 RooRealVar PionDBinaryKaonPionIDSVDonly("PionDBinaryKaonPionIDSVDonly", "", -1.e8, 1.e8);
733 RooRealVar PionDBinaryKaonPionIDnoSVD("PionDBinaryKaonPionIDnoSVD", "", -1.e8, 1.e8);
734
735 RooRealVar PionDBinaryProtonPionIDALL("PionDBinaryProtonPionIDALL", "", -1.e8, 1.e8);
736 RooRealVar PionDBinaryProtonPionIDSVDonly("PionDBinaryProtonPionIDSVDonly", "", -1.e8, 1.e8);
737 RooRealVar PionDBinaryProtonPionIDnoSVD("PionDBinaryProtonPionIDnoSVD", "", -1.e8, 1.e8);
738
739 RooRealVar PionDBinaryElectronPionIDALL("PionDBinaryElectronPionIDALL", "", -1.e8, 1.e8);
740 RooRealVar PionDBinaryElectronPionIDSVDonly("PionDBinaryElectronPionIDSVDonly", "", -1.e8, 1.e8);
741 RooRealVar PionDBinaryElectronPionIDnoSVD("PionDBinaryElectronPionIDnoSVD", "", -1.e8, 1.e8);
742
743 RooRealVar SlowPionKaonIDALL("SlowPionKaonIDALL", "", -1.e8, 1.e8);
744 RooRealVar SlowPionKaonIDSVDonly("SlowPionKaonIDSVDonly", "", -1.e8, 1.e8);
745 RooRealVar SlowPionKaonIDnoSVD("SlowPionKaonIDnoSVD", "", -1.e8, 1.e8);
746
747 RooRealVar SlowPionPionIDALL("SlowPionPionIDALL", "", -1.e8, 1.e8);
748 RooRealVar SlowPionPionIDSVDonly("SlowPionPionIDSVDonly", "", -1.e8, 1.e8);
749 RooRealVar SlowPionPionIDnoSVD("SlowPionPionIDnoSVD", "", -1.e8, 1.e8);
750
751 RooRealVar SlowPionElectronIDALL("SlowPionElectronIDALL", "", -1.e8, 1.e8);
752 RooRealVar SlowPionElectronIDSVDonly("SlowPionElectronIDSVDonly", "", -1.e8, 1.e8);
753 RooRealVar SlowPionElectronIDnoSVD("SlowPionElectronIDnoSVD", "", -1.e8, 1.e8);
754
755 RooRealVar SlowPionProtonIDALL("SlowPionProtonIDALL", "", -1.e8, 1.e8);
756 RooRealVar SlowPionProtonIDSVDonly("SlowPionProtonIDSVDonly", "", -1.e8, 1.e8);
757 RooRealVar SlowPionProtonIDnoSVD("SlowPionProtonIDnoSVD", "", -1.e8, 1.e8);
758
759 RooRealVar SlowPionBinaryPionKaonIDALL("SlowPionBinaryPionKaonIDALL", "", -1.e8, 1.e8);
760 RooRealVar SlowPionBinaryPionKaonIDSVDonly("SlowPionBinaryPionKaonIDSVDonly", "", -1.e8, 1.e8);
761 RooRealVar SlowPionBinaryPionKaonIDnoSVD("SlowPionBinaryPionKaonIDnoSVD", "", -1.e8, 1.e8);
762
763 RooRealVar SlowPionBinaryKaonPionIDALL("SlowPionBinaryKaonPionIDALL", "", -1.e8, 1.e8);
764 RooRealVar SlowPionBinaryKaonPionIDSVDonly("SlowPionBinaryKaonPionIDSVDonly", "", -1.e8, 1.e8);
765 RooRealVar SlowPionBinaryKaonPionIDnoSVD("SlowPionBinaryKaonPionIDnoSVD", "", -1.e8, 1.e8);
766
767 RooRealVar SlowPionBinaryProtonPionIDALL("SlowPionBinaryProtonPionIDALL", "", -1.e8, 1.e8);
768 RooRealVar SlowPionBinaryProtonPionIDSVDonly("SlowPionBinaryProtonPionIDSVDonly", "", -1.e8, 1.e8);
769 RooRealVar SlowPionBinaryProtonPionIDnoSVD("SlowPionBinaryProtonPionIDnoSVD", "", -1.e8, 1.e8);
770
771 RooRealVar SlowPionBinaryElectronPionIDALL("SlowPionBinaryElectronPionIDALL", "", -1.e8, 1.e8);
772 RooRealVar SlowPionBinaryElectronPionIDSVDonly("SlowPionBinaryElectronPionIDSVDonly", "", -1.e8, 1.e8);
773 RooRealVar SlowPionBinaryElectronPionIDnoSVD("SlowPionBinaryElectronPionIDnoSVD", "", -1.e8, 1.e8);
774
775 auto variables = new RooArgSet();
776 variables->add(deltaM);
777 variables->add(KaonMomentum);
778 variables->add(KaonSVDdEdx);
779 variables->add(PionDMomentum);
780 variables->add(PionDSVDdEdx);
781 variables->add(SlowPionMomentum);
782 variables->add(SlowPionSVDdEdx);
783 variables->add(exp);
784 variables->add(run);
785
786 variables->add(KaonKaonIDALL);
787 variables->add(KaonKaonIDSVDonly);
788 variables->add(KaonKaonIDnoSVD);
789 variables->add(KaonPionIDALL);
790 variables->add(KaonPionIDSVDonly);
791 variables->add(KaonPionIDnoSVD);
792 variables->add(KaonProtonIDALL);
793 variables->add(KaonProtonIDSVDonly);
794 variables->add(KaonProtonIDnoSVD);
795 variables->add(KaonElectronIDALL);
796 variables->add(KaonElectronIDSVDonly);
797 variables->add(KaonElectronIDnoSVD);
798 variables->add(KaonBinaryKaonPionIDALL);
799 variables->add(KaonBinaryKaonPionIDSVDonly);
800 variables->add(KaonBinaryKaonPionIDnoSVD);
801 variables->add(KaonBinaryPionKaonIDALL);
802 variables->add(KaonBinaryPionKaonIDSVDonly);
803 variables->add(KaonBinaryPionKaonIDnoSVD);
804 variables->add(KaonBinaryProtonKaonIDALL);
805 variables->add(KaonBinaryProtonKaonIDSVDonly);
806 variables->add(KaonBinaryProtonKaonIDnoSVD);
807 variables->add(KaonBinaryElectronKaonIDALL);
808 variables->add(KaonBinaryElectronKaonIDSVDonly);
809 variables->add(KaonBinaryElectronKaonIDnoSVD);
810
811 variables->add(PionDPionIDALL);
812 variables->add(PionDPionIDSVDonly);
813 variables->add(PionDPionIDnoSVD);
814 variables->add(PionDKaonIDALL);
815 variables->add(PionDKaonIDSVDonly);
816 variables->add(PionDKaonIDnoSVD);
817 variables->add(PionDElectronIDALL);
818 variables->add(PionDElectronIDSVDonly);
819 variables->add(PionDElectronIDnoSVD);
820 variables->add(PionDProtonIDALL);
821 variables->add(PionDProtonIDSVDonly);
822 variables->add(PionDProtonIDnoSVD);
823 variables->add(PionDBinaryPionKaonIDALL);
824 variables->add(PionDBinaryPionKaonIDSVDonly);
825 variables->add(PionDBinaryPionKaonIDnoSVD);
826 variables->add(PionDBinaryKaonPionIDALL);
827 variables->add(PionDBinaryKaonPionIDSVDonly);
828 variables->add(PionDBinaryKaonPionIDnoSVD);
829 variables->add(PionDBinaryProtonPionIDALL);
830 variables->add(PionDBinaryProtonPionIDSVDonly);
831 variables->add(PionDBinaryProtonPionIDnoSVD);
832 variables->add(PionDBinaryElectronPionIDALL);
833 variables->add(PionDBinaryElectronPionIDSVDonly);
834 variables->add(PionDBinaryElectronPionIDnoSVD);
835
836 variables->add(SlowPionPionIDALL);
837 variables->add(SlowPionPionIDSVDonly);
838 variables->add(SlowPionPionIDnoSVD);
839 variables->add(SlowPionKaonIDALL);
840 variables->add(SlowPionKaonIDSVDonly);
841 variables->add(SlowPionKaonIDnoSVD);
842 variables->add(SlowPionElectronIDALL);
843 variables->add(SlowPionElectronIDSVDonly);
844 variables->add(SlowPionElectronIDnoSVD);
845 variables->add(SlowPionProtonIDALL);
846 variables->add(SlowPionProtonIDSVDonly);
847 variables->add(SlowPionProtonIDnoSVD);
848 variables->add(SlowPionBinaryPionKaonIDALL);
849 variables->add(SlowPionBinaryPionKaonIDSVDonly);
850 variables->add(SlowPionBinaryPionKaonIDnoSVD);
851 variables->add(SlowPionBinaryKaonPionIDALL);
852 variables->add(SlowPionBinaryKaonPionIDSVDonly);
853 variables->add(SlowPionBinaryKaonPionIDnoSVD);
854 variables->add(SlowPionBinaryProtonPionIDALL);
855 variables->add(SlowPionBinaryProtonPionIDSVDonly);
856 variables->add(SlowPionBinaryProtonPionIDnoSVD);
857 variables->add(SlowPionBinaryElectronPionIDALL);
858 variables->add(SlowPionBinaryElectronPionIDSVDonly);
859 variables->add(SlowPionBinaryElectronPionIDnoSVD);
860
861 RooDataSet* DstarDataset = new RooDataSet("DstarDataset", "DstarDataset", preselTree.get(), *variables);
862
863 if (DstarDataset->sumEntries() == 0) {
864 B2FATAL("The Dstar dataset is empty, stopping here");
865 }
866
867 RooPlot* DstarFitFrame = DstarDataset->plotOn(deltaM.frame());
868
869 RooRealVar GaussMean("GaussMean", "GaussMean", 0.145, 0.140, 0.150);
870 RooRealVar GaussSigma1("GaussSigma1", "GaussSigma1", 0.01, 1.e-4, 1.0);
871 RooGaussian DstarGauss1("DstarGauss1", "DstarGauss1", deltaM, GaussMean, GaussSigma1);
872 RooRealVar GaussSigma2("GaussSigma2", "GaussSigma2", 0.001, 1.e-4, 1.0);
873 RooGaussian DstarGauss2("DstarGauss2", "DstarGauss2", deltaM, GaussMean, GaussSigma2);
874 RooRealVar fracGaussYield("fracGaussYield", "Fraction of two Gaussians", 0.75, 0.0, 1.0);
875 RooAddPdf DstarSignalPDF("DstarSignalPDF", "DstarGauss1+DstarGauss2", RooArgList(DstarGauss1, DstarGauss2), fracGaussYield);
876
877 RooRealVar dm0Bkg("dm0Bkg", "dm0", 0.13957018, 0.130, 0.140);
878 RooRealVar aBkg("aBkg", "a", -0.0784, -0.08, 3.0);
879 RooRealVar bBkg("bBkg", "b", -0.444713, -0.5, 0.4);
880 RooRealVar cBkg("cBkg", "c", 0.3);
881 RooDstD0BG DstarBkgPDF("DstarBkgPDF", "DstarBkgPDF", deltaM, dm0Bkg, cBkg, aBkg, bBkg);
882 RooRealVar nSignalDstar("nSignalDstar", "signal yield", 0.5 * preselTree->GetEntries(), 0, preselTree->GetEntries());
883 RooRealVar nBkgDstar("nBkgDstar", "background yield", 0.5 * preselTree->GetEntries(), 0, preselTree->GetEntries());
884 RooAddPdf totalPDFDstar("totalPDFDstar", "totalPDFDstar pdf", RooArgList(DstarSignalPDF, DstarBkgPDF),
885 RooArgList(nSignalDstar, nBkgDstar));
886
887 B2INFO("Dstar: Start fitting...");
888 RooFitResult* DstarFitResult = totalPDFDstar.fitTo(*DstarDataset, Save(kTRUE), PrintLevel(-1));
889
890 int status = DstarFitResult->status();
891 int covqual = DstarFitResult->covQual();
892 double diff = nSignalDstar.getValV() + nBkgDstar.getValV() - DstarDataset->sumEntries();
893
894 B2INFO("Dstar: Fit status: " << status << "; covariance quality: " << covqual);
895 // if the fit is not healthy, try again once before giving up, with a slightly different setup:
896 if ((status > 0) || (TMath::Abs(diff) > 1.) || (nSignalDstar.getError() < sqrt(nSignalDstar.getValV()))
897 || (nSignalDstar.getError() > (nSignalDstar.getValV()))) {
898
899 DstarFitResult = totalPDFDstar.fitTo(*DstarDataset, Save(), Strategy(2), Offset(1));
900 status = DstarFitResult->status();
901 covqual = DstarFitResult->covQual();
902 diff = nSignalDstar.getValV() + nBkgDstar.getValV() - DstarDataset->sumEntries();
903 }
904
905 if ((status > 0) || (TMath::Abs(diff) > 1.) || (nSignalDstar.getError() < sqrt(nSignalDstar.getValV()))
906 || (nSignalDstar.getError() > (nSignalDstar.getValV()))) {
907 B2WARNING("Dstar: Fit problem: fit status " << status << "; sum of component yields minus the dataset yield is " << diff <<
908 "; signal yield is " << nSignalDstar.getValV() << ", while its uncertainty is " << nSignalDstar.getError());
909 }
910 if (covqual < 2) {
911 B2INFO("Dstar: Fit warning: covariance quality " << covqual);
912 }
913
914 totalPDFDstar.plotOn(DstarFitFrame, LineColor(TColor::GetColor("#4575b4")));
915
916 double chisquare = DstarFitFrame->chiSquare();
917 B2INFO("Dstar: Fit chi2 = " << chisquare);
918 totalPDFDstar.paramOn(DstarFitFrame, Layout(0.63, 0.96, 0.93), Format("NEU", AutoPrecision(2)));
919 DstarFitFrame->getAttText()->SetTextSize(0.03);
920
921 totalPDFDstar.plotOn(DstarFitFrame, Components("DstarSignalPDF"), LineColor(TColor::GetColor("#d73027")));
922 totalPDFDstar.plotOn(DstarFitFrame, Components("DstarBkgPDF"), LineColor(TColor::GetColor("#fc8d59")));
923 totalPDFDstar.plotOn(DstarFitFrame, LineColor(TColor::GetColor("#4575b4")));
924
925 DstarFitFrame->GetXaxis()->SetTitle("#Deltam [GeV/c^{2}]");
926 TCanvas* canvDstar = new TCanvas("canvDstar", "canvDstar");
927 canvDstar->cd();
928
929 DstarFitFrame->Draw();
930
931 if (m_isMakePlots) {
932 canvDstar->Print("SVDdEdxValidationFitDstar.pdf");
933 TFile DstarFitPlotFile("SVDdEdxValidationDstarFitPlotFile.root", "RECREATE");
934 canvDstar->Write();
935 DstarFitPlotFile.Close();
936 }
937
939
940 RooStats::SPlot* sPlotDatasetDstar = new RooStats::SPlot("sData", "An SPlot", *DstarDataset, &totalPDFDstar,
941 RooArgList(nSignalDstar, nBkgDstar));
942
943 for (int iEvt = 0; iEvt < 5; iEvt++) {
944 if (TMath::Abs(sPlotDatasetDstar->GetSWeight(iEvt, "nSignalDstar") + sPlotDatasetDstar->GetSWeight(iEvt, "nBkgDstar") - 1) > 5.e-3)
945 B2FATAL("Dstar: sPlot error: sum of weights not equal to 1");
946 }
947
948 RooDataSet* DstarDatasetSWeighted = new RooDataSet(DstarDataset->GetName(), DstarDataset->GetTitle(), DstarDataset,
949 *DstarDataset->get());
950
951 RooDataSet::setDefaultStorageType(RooAbsData::Tree);
952 ((RooTreeDataStore*)(DstarDatasetSWeighted->store())->tree())->SetName("treeDstar_sw");
953 TTree* treeDstar_sw = DstarDatasetSWeighted->GetClonedTree();
954
955 B2INFO("Dstar: sPlot done. ");
956
957 return treeDstar_sw;
958}
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:98
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:308
Belle2::SVDdEdxValidationAlgorithm::m_NumEffBins
unsigned int m_NumEffBins
number of bins for the efficiency/fake rate plot
Definition: SVDdEdxValidationAlgorithm.h:100
Belle2::SVDdEdxValidationAlgorithm::m_isMakePlots
bool m_isMakePlots
produce plots for monitoring of the fit quality
Definition: SVDdEdxValidationAlgorithm.h:79
Belle2::SVDdEdxValidationAlgorithm::m_MinEvtsPerTree
int m_MinEvtsPerTree
number of events in TTree below which we don't try to fit
Definition: SVDdEdxValidationAlgorithm.h:96
Belle2::SVDdEdxValidationAlgorithm::LambdaMassFit
TTree * LambdaMassFit(std::shared_ptr< TTree > preselTree)
produce histograms for protons
Definition: SVDdEdxValidationAlgorithm.cc:469
Belle2::SVDdEdxValidationAlgorithm::m_MomHighROC
double m_MomHighROC
upper edge of the momentum interval considered for the ROC curve
Definition: SVDdEdxValidationAlgorithm.h:99
Belle2::SVDdEdxValidationAlgorithm::m_MomHighEff
double m_MomHighEff
upper edge of the momentum interval for the efficiency/fake rate plot
Definition: SVDdEdxValidationAlgorithm.h:101
Belle2::SVDdEdxValidationAlgorithm::m_NumROCpoints
unsigned int m_NumROCpoints
number of points for the ROC curve plotting
Definition: SVDdEdxValidationAlgorithm.h:97
Belle2::SVDdEdxValidationAlgorithm::DstarMassFit
TTree * DstarMassFit(std::shared_ptr< TTree > preselTree)
produce histograms for K/pi(/mu)
Definition: SVDdEdxValidationAlgorithm.cc:661
Belle2::SVDdEdxValidationAlgorithm::calibrate
virtual EResult calibrate() override
run algorithm on data
Definition: SVDdEdxValidationAlgorithm.cc:51
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:138
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