Belle II Software development
SVDdEdxCalibrationAlgorithm.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/SVDdEdxCalibrationAlgorithm.h>
10#include <svd/dbobjects/SVDdEdxPDFs.h>
11
12#include <TROOT.h>
13#include <TStyle.h>
14#include <TMath.h>
15#include <TFile.h>
16#include <TTree.h>
17#include <TColor.h>
18#include <TLegend.h>
19#include <TCanvas.h>
20#include <TH1D.h>
21#include <TAxis.h>
22#include <TFitResult.h>
23#include <TKDTreeBinning.h>
24
25#include <RooDataSet.h>
26#include <RooRealVar.h>
27#include <RooAddPdf.h>
28#include <RooGaussian.h>
29#include <RooChebychev.h>
30#include <RooBifurGauss.h>
31#include <RooDstD0BG.h>
32#include <RooAbsDataStore.h>
33#include <RooTreeDataStore.h>
34#include <RooMsgService.h>
35#include <RooStats/SPlot.h>
36#include <Math/MinimizerOptions.h>
37
38using namespace RooFit;
39using namespace Belle2;
40
41SVDdEdxCalibrationAlgorithm::SVDdEdxCalibrationAlgorithm() : CalibrationAlgorithm("SVDdEdxCollector"),
42 m_isMakePlots(true)
43{
44 setDescription("SVD dE/dx calibration algorithm");
45}
46
47/* Main calibration method */
48CalibrationAlgorithm::EResult SVDdEdxCalibrationAlgorithm::calibrate()
49{
50 gROOT->SetBatch(true);
51
52 const auto exprun = getRunList()[0];
53 B2INFO("ExpRun used for calibration: " << exprun.first << " " << exprun.second);
54
55 auto payload = new Belle2::SVDdEdxPDFs();
56
57 // Get data objects
58 auto ttreeLambda = getObjectPtr<TTree>("Lambda");
59 auto ttreeDstar = getObjectPtr<TTree>("Dstar");
60 auto ttreeGamma = getObjectPtr<TTree>("Gamma");
61 auto ttreeGeneric = getObjectPtr<TTree>("Generic");
62
63 if ((ttreeLambda->GetEntries() < m_MinEvtsPerTree) || (ttreeDstar->GetEntries() < m_MinEvtsPerTree)
64 || (ttreeGamma->GetEntries() < m_MinEvtsPerTree)) {
65 B2WARNING("Not enough data for calibration.");
66 return c_NotEnoughData;
67 }
68
69 // call the calibration function
70 std::unique_ptr<TList> GeneratedList = GenerateNewHistograms(ttreeLambda, ttreeDstar, ttreeGamma, ttreeGeneric);
71
72 TH2F* histoE = static_cast<TH2F*>(GeneratedList->FindObject("Electron2DHistogramNew"));
73 TH2F* histoMu = static_cast<TH2F*>(GeneratedList->FindObject("Muon2DHistogramNew"));
74 TH2F* histoPi = static_cast<TH2F*>(GeneratedList->FindObject("Pion2DHistogramNew"));
75 TH2F* histoK = static_cast<TH2F*>(GeneratedList->FindObject("Kaon2DHistogramNew"));
76 TH2F* histoP = static_cast<TH2F*>(GeneratedList->FindObject("Proton2DHistogramNew"));
77 TH2F* histoDeut = static_cast<TH2F*>(GeneratedList->FindObject("Deuteron2DHistogramNew"));
78
79 std::vector<double> pbins = CreatePBinningScheme();
80 TH2F hEmpty("hEmpty", "A histogram returned if we cannot calibrate", m_numPBins, pbins.data(), m_numDEdxBins, 0, m_dedxCutoff);
81 for (int pbin = 0; pbin <= m_numPBins + 1; pbin++) {
82 for (int dedxbin = 0; dedxbin <= m_numDEdxBins + 1; dedxbin++) {
83 hEmpty.SetBinContent(pbin, dedxbin, 0.01);
84 };
85 }
86
87 B2INFO("Histograms are ready, proceed to creating the payload object...");
88 std::vector<TH2F*> hDedxPDFs(6);
89
90 std::array<std::string, 6> part = {"Electron", "Muon", "Pion", "Kaon", "Proton", "Deuteron"};
91
92 std::unique_ptr<TCanvas> candEdx(new TCanvas("candEdx", "SVD dEdx payloads", 1200, 700));
93 candEdx->Divide(3, 2);
94 gStyle->SetOptStat(11);
95
96 for (bool trunmean : {false, true}) {
97 for (int iPart = 0; iPart < 6; iPart++) {
98 if (iPart == 0 && trunmean) {
99 hDedxPDFs[iPart] = histoE;
100 hDedxPDFs[iPart]->SetName("hist_d1_11_trunc");
101 } else if (iPart == 1 && trunmean) {
102 hDedxPDFs[iPart] = histoMu;
103 hDedxPDFs[iPart]->SetName("hist_d1_13_trunc");
104 } else if (iPart == 2 && trunmean) {
105 hDedxPDFs[iPart] = histoPi;
106 hDedxPDFs[iPart]->SetName("hist_d1_211_trunc");
107 } else if (iPart == 3 && trunmean) {
108 hDedxPDFs[iPart] = histoK;
109 hDedxPDFs[iPart]->SetName("hist_d1_321_trunc");
110 } else if (iPart == 4 && trunmean) {
111 hDedxPDFs[iPart] = histoP;
112 hDedxPDFs[iPart]->SetName("hist_d1_2212_trunc");
113 } else if (iPart == 5 && trunmean) {
114 hDedxPDFs[iPart] = histoDeut;
115 hDedxPDFs[iPart]->SetName("hist_d1_1000010020_trunc");
116 } else if (iPart == 0 && !trunmean) {
117 hDedxPDFs[iPart] = &hEmpty;
118 hDedxPDFs[iPart]->SetName("hist_d1_11");
119 } else if (iPart == 1 && !trunmean) {
120 hDedxPDFs[iPart] = &hEmpty;
121 hDedxPDFs[iPart]->SetName("hist_d1_13");
122 } else if (iPart == 2 && !trunmean) {
123 hDedxPDFs[iPart] = &hEmpty;
124 hDedxPDFs[iPart]->SetName("hist_d1_211");
125 } else if (iPart == 3 && !trunmean) {
126 hDedxPDFs[iPart] = &hEmpty;
127 hDedxPDFs[iPart]->SetName("hist_d1_321");
128 } else if (iPart == 4 && !trunmean) {
129 hDedxPDFs[iPart] = &hEmpty;
130 hDedxPDFs[iPart]->SetName("hist_d1_2212");
131 } else if (iPart == 5 && !trunmean) {
132 hDedxPDFs[iPart] = &hEmpty;
133 hDedxPDFs[iPart]->SetName("hist_d1_1000010020");
134 } else
135 hDedxPDFs[iPart] = &hEmpty;
136 payload->setPDF(*hDedxPDFs[iPart], iPart, trunmean);
137
138 candEdx->cd(iPart + 1);
139 hDedxPDFs[iPart]->SetTitle(Form("%s; p(GeV/c) of %s; dE/dx", hDedxPDFs[iPart]->GetTitle(), part[iPart].data()));
140 hDedxPDFs[iPart]->DrawCopy("colz");
141 }
142
143 if (m_isMakePlots) {
144 candEdx->SaveAs("PlotsSVDdEdxPDFs_wTruncMean.pdf");
145 std::unique_ptr<TList> l(new TList());
146 for (int iPart = 0; iPart < 6; iPart++) {
147 l->Add(hDedxPDFs[iPart]);
148 }
149
150 TFile SVDdEdxPDFsPlotFile("PlotsSVDdEdxPDFs_wTruncMean.root", "RECREATE");
151 l->Write("histlist", TObject::kSingleKey);
152 SVDdEdxPDFsPlotFile.Close();
153 }
154
155 // candEdx->SetTitle(Form("Likehood dist. of charged particles from %s, trunmean = %s", idet.data(), check.str().data()));
156 }
157
158 saveCalibration(payload, "SVDdEdxPDFs");
159 B2INFO("SVD dE/dx calibration done!");
160
161 return c_OK;
162}
163
164TTree* SVDdEdxCalibrationAlgorithm::LambdaMassFit(std::shared_ptr<TTree> preselTree)
165{
166 B2INFO("Configuring the Lambda fit...");
167 gROOT->SetBatch(true);
168 RooMsgService::instance().setGlobalKillBelow(RooFit::WARNING);
169
170 RooRealVar InvM("InvM", "m(p^{+}#pi^{-})", 1.1, 1.13, "GeV/c^{2}");
171
172 RooRealVar ProtonMomentum("ProtonMomentum", "momentum for p", -1.e8, 1.e8);
173 RooRealVar ProtonSVDdEdxTrackMomentum("ProtonSVDdEdxTrackMomentum", "momentum for p", -1.e8, 1.e8);
174 RooRealVar ProtonSVDdEdx("ProtonSVDdEdx", "", -1.e8, 1.e8);
175 RooRealVar ProtonSVDdEdxTrackCosTheta("ProtonSVDdEdxTrackCosTheta", "", -10., 10.);
176
177 RooRealVar exp("exp", "experiment number", 0, 1.e5);
178 RooRealVar run("run", "run number", 0, 1.e7);
179
180 auto variables = new RooArgSet();
181
182 variables->add(InvM);
183
184 variables->add(ProtonMomentum);
185 variables->add(ProtonSVDdEdxTrackMomentum);
186 variables->add(ProtonSVDdEdx);
187 variables->add(ProtonSVDdEdxTrackCosTheta);
188 variables->add(exp);
189 variables->add(run);
190
191 RooDataSet* LambdaDataset = new RooDataSet("LambdaDataset", "LambdaDataset", *variables, Import(*preselTree));
192
193 if (LambdaDataset->sumEntries() == 0) {
194 B2FATAL("The Lambda dataset is empty, stopping here");
195 }
196
197 // the signal PDF; might be revisited at a later point
198
199 RooRealVar GaussMean("GaussMean", " GaussMean", 1.116, 1.111, 1.12);
200 RooRealVar GaussSigma("GaussSigma", "#sigma_{1}", 3.e-3, 3.e-5, 10.e-3);
201 RooGaussian LambdaGauss("LambdaGauss", "LambdaGauss", InvM, GaussMean, GaussSigma);
202
203 /* temporary RooRealVar sigmaBifurGaussL1 and sigmaBifurGaussR1 to replace
204 * RooRealVar resolutionParamL("resolutionParamL", "resolutionParamL", 0.4, 5.e-4, 1.0);
205 * RooRealVar resolutionParamR("resolutionParamR", "resolutionParamR", 0.4, 5.e-4, 1.0);
206 * RooFormulaVar sigmaBifurGaussL1("sigmaBifurGaussL1", "resolutionParamL*GaussSigma", RooArgSet(resolutionParamL, GaussSigma));
207 * RooFormulaVar sigmaBifurGaussR1("sigmaBifurGaussR1", "resolutionParamR*GaussSigma", RooArgSet(resolutionParamR, GaussSigma));
208 */
209 RooRealVar sigmaBifurGaussL1("sigmaBifurGaussL1", "sigma left", 0.4 * 3.e-3, 3.e-5, 10.e-3);
210 RooRealVar sigmaBifurGaussR1("sigmaBifurGaussR1", "sigma right", 0.4 * 3.e-3, 3.e-5, 10.e-3);
211 RooBifurGauss LambdaBifurGauss("LambdaBifurGauss", "LambdaBifurGauss", InvM, GaussMean, sigmaBifurGaussL1, sigmaBifurGaussR1);
212
213 /* temporary RooRealVar sigmaBifurGaussL2 to replace
214 * RooRealVar resolutionParam2("resolutionParam2", "resolutionParam2", 0.2, 5.e-4, 1.0);
215 * sigmaBifurGaussL2("sigmaBifurGaussL2", "resolutionParam2*GaussSigma", RooArgSet(resolutionParam2, GaussSigma));
216 */
217 RooRealVar sigmaBifurGaussL2("sigmaBifurGaussL2", "sigmaBifurGaussL2", 0.2 * 3.e-3, 3.e-5, 10.e-3);
218 RooGaussian LambdaBifurGauss2("LambdaBifurGauss2", "LambdaBifurGauss2", InvM, GaussMean, sigmaBifurGaussL2);
219
220 RooRealVar fracBifurGaussYield("fracBifurGaussYield", "fracBifurGaussYield", 0.3, 5.e-4, 1.0);
221 RooRealVar fracGaussYield("fracGaussYield", "fracGaussYield", 0.8, 5.e-4, 1.0);
222
223 RooAddPdf LambdaCombinedBifurGauss("LambdaCombinedBifurGauss", "LambdaBifurGauss + LambdaBifurGauss2 ", RooArgList(LambdaBifurGauss,
224 LambdaBifurGauss2), RooArgList(fracBifurGaussYield));
225
226 RooAddPdf LambdaSignalPDF("LambdaSignalPDF", "LambdaCombinedBifurGauss + LambdaGauss", RooArgList(LambdaCombinedBifurGauss,
227 LambdaGauss), RooArgList(fracGaussYield));
228
229 // Background PDF
230 RooRealVar BkgPolyCoef0("BkgPolyCoef0", "BkgPolyCoef0", 0.1, 0., 1.5);
231 RooRealVar BkgPolyCoef1("BkgPolyCoef1", "BkgPolyCoef1", -0.5, -1.5, -1.e-3);
232 RooChebychev BkgPolyPDF("BkgPolyPDF", "BkgPolyPDF", InvM, RooArgList(BkgPolyCoef0, BkgPolyCoef1));
233
234 RooRealVar nSignalLambda("nSignalLambda", "nSignalLambda", 0.6 * preselTree->GetEntries(), 0., 0.99 * preselTree->GetEntries());
235 RooRealVar nBkgLambda("nBkgLambda", "nBkgLambda", 0.4 * preselTree->GetEntries(), 0., 0.99 * preselTree->GetEntries());
236 RooAddPdf totalPDFLambda("totalPDFLambda", "totalPDFLambda pdf", RooArgList(LambdaSignalPDF, BkgPolyPDF),
237 RooArgList(nSignalLambda, nBkgLambda));
238
239 B2INFO("Lambda: Start fitting...");
240 RooFitResult* LambdaFitResult = totalPDFLambda.fitTo(*LambdaDataset, Save(kTRUE), PrintLevel(-1));
241
242 int status = LambdaFitResult->status();
243 int covqual = LambdaFitResult->covQual();
244 double diff = nSignalLambda.getValV() + nBkgLambda.getValV() - LambdaDataset->sumEntries();
245
246 B2INFO("Lambda: Fit status: " << status << "; covariance quality: " << covqual);
247 // if the fit is not healthy, try again once before giving up, with a slightly different setup:
248 if ((status > 0) || (TMath::Abs(diff) > 1.) || (nSignalLambda.getError() < sqrt(nSignalLambda.getValV()))
249 || (nSignalLambda.getError() > (nSignalLambda.getValV()))) {
250
251 LambdaFitResult = totalPDFLambda.fitTo(*LambdaDataset, Save(), Strategy(2), Offset(1));
252 status = LambdaFitResult->status();
253 covqual = LambdaFitResult->covQual();
254 diff = nSignalLambda.getValV() + nBkgLambda.getValV() - LambdaDataset->sumEntries();
255 B2INFO("Lambda: updated fit status: " << status << "; covariance quality: " << covqual);
256 }
257
258 if ((status > 0) || (TMath::Abs(diff) > 1.) || (nSignalLambda.getError() < sqrt(nSignalLambda.getValV()))
259 || (nSignalLambda.getError() > (nSignalLambda.getValV()))) {
260 B2WARNING("Lambda: Fit problem: fit status " << status << "; sum of component yields minus the dataset yield is " << diff <<
261 "; signal yield is " << nSignalLambda.getValV() << ", while its uncertainty is " << nSignalLambda.getError());
262 }
263 if (covqual < 2) {
264 B2INFO("Lambda: Fit warning: covariance quality " << covqual);
265 }
266
267 if (m_isMakePlots) {
268 std::unique_ptr<TCanvas> canvLambda(new TCanvas("canvLambda", "canvLambda"));
269 canvLambda->cd();
270 RooPlot* LambdaFitFrame = LambdaDataset->plotOn(InvM.frame(130));
271 totalPDFLambda.plotOn(LambdaFitFrame, LineColor(TColor::GetColor("#4575b4")));
272
273 double chisquare = LambdaFitFrame->chiSquare();
274 B2INFO("Lambda: Fit chi2 = " << chisquare);
275 totalPDFLambda.paramOn(LambdaFitFrame, Layout(0.6, 0.96, 0.93), Format("NEU", AutoPrecision(2)));
276 LambdaFitFrame->getAttText()->SetTextSize(0.03);
277
278 totalPDFLambda.plotOn(LambdaFitFrame, Components("LambdaSignalPDF"), LineColor(TColor::GetColor("#d73027")));
279 totalPDFLambda.plotOn(LambdaFitFrame, Components("BkgPolyPDF"), LineColor(TColor::GetColor("#fc8d59")));
280 totalPDFLambda.plotOn(LambdaFitFrame, LineColor(TColor::GetColor("#4575b4")));
281
282 LambdaFitFrame->GetXaxis()->SetTitle("m(p#pi^{-}) (GeV/c^{2})");
283
284 LambdaFitFrame->Draw();
285
286
287 canvLambda->Print("SVDdEdxCalibrationFitLambda.pdf");
288 TFile LambdaFitPlotFile("SVDdEdxCalibrationLambdaFitPlotFile.root", "RECREATE");
289 canvLambda->Write();
290 LambdaFitPlotFile.Close();
291 }
292 RooStats::SPlot* sPlotDatasetLambda = new RooStats::SPlot("sData", "An SPlot", *LambdaDataset, &totalPDFLambda,
293 RooArgList(nSignalLambda, nBkgLambda));
294
295 for (int iEvt = 0; iEvt < 5; iEvt++) {
296 if (TMath::Abs(sPlotDatasetLambda->GetSWeight(iEvt, "nSignalLambda") + sPlotDatasetLambda->GetSWeight(iEvt,
297 "nBkgLambda") - 1) > 5.e-3)
298 B2FATAL("Lambda: sPlot error: sum of weights not equal to 1");
299 }
300
301 TTree* treeLambdaSWeighted = LambdaDataset->GetClonedTree();
302 treeLambdaSWeighted->SetName("treeLambdaSWeighted");
303
304 B2INFO("Lambda: sPlot done. Proceed to histogramming");
305 return treeLambdaSWeighted;
306}
307
308std::unique_ptr<TList> SVDdEdxCalibrationAlgorithm::LambdaHistogramming(TTree* inputTree)
309{
310 gROOT->SetBatch(true);
311 inputTree->SetEstimate(-1);
312 std::vector<double> pbins = CreatePBinningScheme();
313
314 TH2F* hLambdaPMomentum = new TH2F("hist_d1_2212_truncMomentum", "hist_d1_2212_trunc;Momentum [GeV/c];dEdx [arb. units]",
315 m_numPBins, pbins.data(), m_numDEdxBins, 0,
316 m_dedxCutoff);
317
318 inputTree->Draw("ProtonSVDdEdx:ProtonSVDdEdxTrackMomentum>>hist_d1_2212_truncMomentum",
319 "nSignalLambda_sw * (ProtonSVDdEdx>0) * (ProtonSVDdEdxTrackMomentum>0.13)", "goff");
320
321// create isopopulated beta*gamma binning
322 inputTree->Draw(Form("ProtonSVDdEdxTrackMomentum/%f", m_ProtonPDGMass), "", "goff",
323 ((inputTree->GetEntries()) / m_numBGBins)*m_numBGBins);
324 double* ProtonMomentumDataset = inputTree->GetV1();
325
326 TKDTreeBinning* kdBinsP = new TKDTreeBinning(((inputTree->GetEntries()) / m_numBGBins)*m_numBGBins, 1, ProtonMomentumDataset,
327 m_numBGBins);
328 const double* binsMinEdgesP_pointer = kdBinsP->SortOneDimBinEdges();
329 double* binsMinEdgesP = const_cast<double*>(binsMinEdgesP_pointer);
330
331
332 binsMinEdgesP[0] = 0.1;
333 binsMinEdgesP[m_numBGBins + 1] = 50.;
334
335
336 TH2F* hLambdaPBetaGamma = new TH2F("hist_d1_2212_truncBetaGamma", "hist_d1_2212_truncBetaGamma;#beta*#gamma;dEdx [arb. units]",
337 m_numBGBins, binsMinEdgesP, m_numDEdxBins,
338 0,
339 m_dedxMaxPossible);
340
341 inputTree->Draw(Form("ProtonSVDdEdx:ProtonSVDdEdxTrackMomentum/%f>>hist_d1_2212_truncBetaGamma", m_ProtonPDGMass),
342 "nSignalLambda_sw * (ProtonSVDdEdx>0) * (ProtonSVDdEdxTrackMomentum>0.13) * (ProtonSVDdEdx>1.2e6 - 1.e6*ProtonSVDdEdxTrackMomentum)",
343 "goff");
344
345 // produce the 1D profile
346 // momentum: for data-MC comparisons
347
348 TH1D* ProtonProfileMomentum = static_cast<TH1D*>(hLambdaPMomentum->ProfileX("ProtonProfileMomentum"));
349 ProtonProfileMomentum->SetTitle("ProtonProfile");
350 ProtonProfileMomentum->GetYaxis()->SetRangeUser(0, m_dedxCutoff);
351 ProtonProfileMomentum->GetXaxis()->SetTitle("Momentum, GeV/c");
352 ProtonProfileMomentum->GetYaxis()->SetTitle("dE/dx");
353 ProtonProfileMomentum->SetLineColor(kRed);
354
355 //beta*gamma: for the fit
356 TH1D* ProtonProfileBetaGamma = static_cast<TH1D*>(hLambdaPBetaGamma->ProfileX("ProtonProfileBetaGamma"));
357 if (m_CustomProfile) {
358 ProtonProfileBetaGamma = PrepareProfile(hLambdaPBetaGamma, "ProtonProfileBetaGamma");
359 }
360 ProtonProfileBetaGamma->SetTitle("ProtonProfile");
361 ProtonProfileBetaGamma->GetYaxis()->SetRangeUser(0, m_dedxCutoff);
362 ProtonProfileBetaGamma->GetXaxis()->SetTitle("#beta*#gamma");
363 ProtonProfileBetaGamma->GetYaxis()->SetTitle("dE/dx");
364 ProtonProfileBetaGamma->SetLineColor(kRed);
365
366
367 // for each momentum bin, normalize the pdf
368 hLambdaPMomentum = Normalise2DHisto(hLambdaPMomentum);
369
370 std::unique_ptr<TList> histList(new TList);
371 histList->Add(ProtonProfileMomentum);
372 histList->Add(ProtonProfileBetaGamma);
373 histList->Add(hLambdaPMomentum);
374
375 if (m_isMakePlots) {
376 TFile LambdaHistogrammingPlotFile("SVDdEdxCalibrationLambdaHistogramming.root", "RECREATE");
377 histList->Write();
378 LambdaHistogrammingPlotFile.Close();
379 }
380
381 return histList;
382}
383
384TTree* SVDdEdxCalibrationAlgorithm::DstarMassFit(std::shared_ptr<TTree> preselTree)
385{
386 B2INFO("Configuring the Dstar fit...");
387 gROOT->SetBatch(true);
388 RooMsgService::instance().setGlobalKillBelow(RooFit::WARNING);
389
390 RooRealVar deltaM("deltaM", "m(D*)-m(D^{0})", 0.139545, 0.151, "GeV/c^{2}");
391
392 RooRealVar KaonMomentum("KaonMomentum", "momentum for Kaon (GeV)", -1.e8, 1.e8);
393 RooRealVar KaonSVDdEdxTrackMomentum("KaonSVDdEdxTrackMomentum", "momentum for Kaon (GeV), from the track", -1.e8, 1.e8);
394 RooRealVar KaonSVDdEdx("KaonSVDdEdx", "", -1.e8, 1.e8);
395 RooRealVar PionDMomentum("PionDMomentum", "momentum for pion (GeV)", -1.e8, 1.e8);
396 RooRealVar PionDSVDdEdxTrackMomentum("PionDSVDdEdxTrackMomentum", "momentum for pion (GeV), from the track", -1.e8, 1.e8);
397 RooRealVar PionDSVDdEdx("PionDSVDdEdx", "", -1.e8, 1.e8);
398 RooRealVar SlowPionMomentum("SlowPionMomentum", "momentum for slow pion (GeV)", -1.e8, 1.e8);
399 RooRealVar SlowPionSVDdEdxTrackMomentum("SlowPionSVDdEdxTrackMomentum", "momentum for slow pion (GeV), from the track", -1.e8,
400 1.e8);
401 RooRealVar SlowPionSVDdEdx("SlowPionSVDdEdx", "", -1.e8, 1.e8);
402
403 RooRealVar exp("exp", "experiment number", 0, 1.e5);
404 RooRealVar run("run", "run number", 0, 1.e8);
405 RooRealVar event("event", "event number", 0, 1.e10);
406
407 auto variables = new RooArgSet();
408 variables->add(deltaM);
409 variables->add(KaonMomentum);
410 variables->add(KaonSVDdEdxTrackMomentum);
411 variables->add(KaonSVDdEdx);
412 variables->add(PionDMomentum);
413 variables->add(PionDSVDdEdxTrackMomentum);
414 variables->add(PionDSVDdEdx);
415 variables->add(SlowPionMomentum);
416 variables->add(SlowPionSVDdEdxTrackMomentum);
417 variables->add(SlowPionSVDdEdx);
418 variables->add(exp);
419 variables->add(run);
420 variables->add(event);
421
422 RooDataSet* DstarDataset = new RooDataSet("DstarDataset", "DstarDataset", *variables, Import(*preselTree));
423
424 if (DstarDataset->sumEntries() == 0) {
425 B2FATAL("The Dstar dataset is empty, stopping here");
426 }
427
428 RooPlot* DstarFitFrame = DstarDataset->plotOn(deltaM.frame());
429
430 RooRealVar GaussMean("GaussMean", "GaussMean", 0.145, 0.140, 0.150);
431 RooRealVar GaussSigma1("GaussSigma1", "GaussSigma1", 0.01, 1.e-4, 1.0);
432 RooGaussian DstarGauss1("DstarGauss1", "DstarGauss1", deltaM, GaussMean, GaussSigma1);
433 RooRealVar GaussSigma2("GaussSigma2", "GaussSigma2", 0.001, 1.e-4, 1.0);
434 RooGaussian DstarGauss2("DstarGauss2", "DstarGauss2", deltaM, GaussMean, GaussSigma2);
435 RooRealVar fracGaussYield("fracGaussYield", "Fraction of two Gaussians", 0.75, 0.0, 1.0);
436 RooAddPdf DstarSignalPDF("DstarSignalPDF", "DstarGauss1+DstarGauss2", RooArgList(DstarGauss1, DstarGauss2), fracGaussYield);
437
438 RooRealVar dm0Bkg("dm0Bkg", "dm0", 0.13957018, 0.130, 0.140);
439 RooRealVar aBkg("aBkg", "a", -0.0784, -0.08, 3.0);
440 RooRealVar bBkg("bBkg", "b", -0.444713, -0.5, 0.4);
441 RooRealVar cBkg("cBkg", "c", 0.3);
442 RooDstD0BG DstarBkgPDF("DstarBkgPDF", "DstarBkgPDF", deltaM, dm0Bkg, cBkg, aBkg, bBkg);
443 RooRealVar nSignalDstar("nSignalDstar", "signal yield", 0.5 * preselTree->GetEntries(), 0, preselTree->GetEntries());
444 RooRealVar nBkgDstar("nBkgDstar", "background yield", 0.5 * preselTree->GetEntries(), 0, preselTree->GetEntries());
445 RooAddPdf totalPDFDstar("totalPDFDstar", "totalPDFDstar pdf", RooArgList(DstarSignalPDF, DstarBkgPDF),
446 RooArgList(nSignalDstar, nBkgDstar));
447
448 B2INFO("Dstar: Start fitting...");
449 RooFitResult* DstarFitResult = totalPDFDstar.fitTo(*DstarDataset, Save(kTRUE), PrintLevel(-1));
450
451 int status = DstarFitResult->status();
452 int covqual = DstarFitResult->covQual();
453 double diff = nSignalDstar.getValV() + nBkgDstar.getValV() - DstarDataset->sumEntries();
454
455 B2INFO("Dstar: Fit status: " << status << "; covariance quality: " << covqual);
456 // if the fit is not healthy, try again once before giving up, with a slightly different setup:
457 if ((status > 0) || (TMath::Abs(diff) > 1.) || (nSignalDstar.getError() < sqrt(nSignalDstar.getValV()))
458 || (nSignalDstar.getError() > (nSignalDstar.getValV()))) {
459
460 DstarFitResult = totalPDFDstar.fitTo(*DstarDataset, Save(), Strategy(2), Offset(1));
461 status = DstarFitResult->status();
462 covqual = DstarFitResult->covQual();
463 diff = nSignalDstar.getValV() + nBkgDstar.getValV() - DstarDataset->sumEntries();
464 B2INFO("Dstar: Updated fit status: " << status << "; covariance quality: " << covqual);
465 }
466
467 if ((status > 0) || (TMath::Abs(diff) > 1.) || (nSignalDstar.getError() < sqrt(nSignalDstar.getValV()))
468 || (nSignalDstar.getError() > (nSignalDstar.getValV()))) {
469 B2WARNING("Dstar: Fit problem: fit status " << status << "; sum of component yields minus the dataset yield is " << diff <<
470 "; signal yield is " << nSignalDstar.getValV() << ", while its uncertainty is " << nSignalDstar.getError());
471 }
472 if (covqual < 2) {
473 B2INFO("Dstar: Fit warning: covariance quality " << covqual);
474 }
475
476 totalPDFDstar.plotOn(DstarFitFrame, LineColor(TColor::GetColor("#4575b4")));
477
478 double chisquare = DstarFitFrame->chiSquare();
479 B2INFO("Dstar: Fit chi2 = " << chisquare);
480 totalPDFDstar.paramOn(DstarFitFrame, Layout(0.63, 0.96, 0.93), Format("NEU", AutoPrecision(2)));
481 DstarFitFrame->getAttText()->SetTextSize(0.03);
482
483 totalPDFDstar.plotOn(DstarFitFrame, Components("DstarSignalPDF"), LineColor(TColor::GetColor("#d73027")));
484 totalPDFDstar.plotOn(DstarFitFrame, Components("DstarBkgPDF"), LineColor(TColor::GetColor("#fc8d59")));
485 totalPDFDstar.plotOn(DstarFitFrame, LineColor(TColor::GetColor("#4575b4")));
486
487 DstarFitFrame->GetXaxis()->SetTitle("#Deltam [GeV/c^{2}]");
488 if (m_isMakePlots) {
489 std::unique_ptr<TCanvas> canvDstar(new TCanvas("canvDstar", "canvDstar"));
490 canvDstar->cd();
491
492 DstarFitFrame->Draw();
493
494 canvDstar->Print("SVDdEdxCalibrationFitDstar.pdf");
495 TFile DstarFitPlotFile("SVDdEdxCalibrationDstarFitPlotFile.root", "RECREATE");
496 canvDstar->Write();
497 DstarFitPlotFile.Close();
498 }
499
501
502 RooStats::SPlot* sPlotDatasetDstar = new RooStats::SPlot("sData", "An SPlot", *DstarDataset, &totalPDFDstar,
503 RooArgList(nSignalDstar, nBkgDstar));
504
505 for (int iEvt = 0; iEvt < 5; iEvt++) {
506 if (TMath::Abs(sPlotDatasetDstar->GetSWeight(iEvt, "nSignalDstar") + sPlotDatasetDstar->GetSWeight(iEvt, "nBkgDstar") - 1) > 5.e-3)
507 B2FATAL("Dstar: sPlot error: sum of weights not equal to 1");
508 }
509
510 TTree* treeDstarSWeighted = DstarDataset->GetClonedTree();
511 treeDstarSWeighted->SetName("treeDstarSWeighted");
512
513 B2INFO("Dstar: sPlot done. Proceed to histogramming");
514 return treeDstarSWeighted;
515}
516
517std::unique_ptr<TList> SVDdEdxCalibrationAlgorithm::DstarHistogramming(TTree* inputTree)
518{
519 gROOT->SetBatch(true);
520 inputTree->SetEstimate(-1);
521 std::vector<double> pbins = CreatePBinningScheme();
522
523 TH2F* hDstarKMomentum = new TH2F("hist_d1_321_truncMomentum", "hist_d1_321_trunc;Momentum [GeV/c];dEdx [arb. units]", m_numPBins,
524 pbins.data(),
525 m_numDEdxBins, 0, m_dedxCutoff);
526 // the pion payload
527 TH2F* hDstarPiMomentum = new TH2F("hist_d1_211_truncMomentum", "hist_d1_211_trunc;Momentum [GeV/c];dEdx [arb. units]", m_numPBins,
528 pbins.data(),
529 m_numDEdxBins, 0, m_dedxCutoff);
530
531 inputTree->Draw("KaonSVDdEdx:KaonSVDdEdxTrackMomentum>>hist_d1_321_truncMomentum", "nSignalDstar_sw * (KaonSVDdEdx>0)", "goff");
532 // the pion one will be built from both pions in the Dstar decay tree
533 TH2F* hDstarPiPart1Momentum = static_cast<TH2F*>(hDstarPiMomentum->Clone("hist_d1_211_truncPart1Momentum"));
534 TH2F* hDstarPiPart2Momentum = static_cast<TH2F*>(hDstarPiMomentum->Clone("hist_d1_211_truncPart2Momentum"));
535
536 inputTree->Draw("PionDSVDdEdx:PionDSVDdEdxTrackMomentum>>hist_d1_211_truncPart1Momentum", "nSignalDstar_sw * (PionDSVDdEdx>0)",
537 "goff");
538 inputTree->Draw("SlowPionSVDdEdx:SlowPionSVDdEdxTrackMomentum>>hist_d1_211_truncPart2Momentum",
539 "nSignalDstar_sw * (SlowPionSVDdEdx>0)",
540 "goff");
541 hDstarPiMomentum->Add(hDstarPiPart1Momentum);
542 hDstarPiMomentum->Add(hDstarPiPart2Momentum);
543
544 inputTree->Draw(Form("KaonSVDdEdxTrackMomentum/%f", m_KaonPDGMass), "", "goff",
545 ((inputTree->GetEntries()) / m_numBGBins)*m_numBGBins);
546 double* KaonMomentumDataset = inputTree->GetV1();
547 TKDTreeBinning* kdBinsK = new TKDTreeBinning(((inputTree->GetEntries()) / m_numBGBins)*m_numBGBins, 1, KaonMomentumDataset,
548 m_numBGBins);
549 const double* binsMinEdgesKOriginal = kdBinsK->SortOneDimBinEdges();
550 double* binsMinEdgesK = const_cast<double*>(binsMinEdgesKOriginal);
551 binsMinEdgesK[0] = 0.1;
552 binsMinEdgesK[m_numBGBins + 1] = 50.;
553
554// get a distribution that contains both pions to get a typical kinematics for the binning scheme
555 inputTree->Draw(Form("SlowPionSVDdEdxTrackMomentum/%f * (event %% 2 == 0) + PionDSVDdEdxTrackMomentum/%f * (event %% 2 ==1)",
556 m_PionPDGMass, m_PionPDGMass), "", "goff",
557 ((inputTree->GetEntries()) / m_numBGBins)*m_numBGBins);
558 double* PionMomentumDataset = inputTree->GetV1();
559
560 TKDTreeBinning* kdBinsPi = new TKDTreeBinning(((inputTree->GetEntries()) / m_numBGBins)*m_numBGBins, 1, PionMomentumDataset,
561 m_numBGBins);
562 const double* binsMinEdgesPiOriginal = kdBinsPi->SortOneDimBinEdges();
563 double* binsMinEdgesPi = const_cast<double*>(binsMinEdgesPiOriginal);
564 binsMinEdgesPi[0] = 0.1;
565 binsMinEdgesPi[m_numBGBins + 1] = 50.;
566
567 TH2F* hDstarKBetaGamma = new TH2F("hist_d1_321_truncBetaGamma", "hist_d1_321_truncBetaGamma;#beta*#gamma;dEdx [arb. units]",
568 m_numBGBins,
569 binsMinEdgesK,
570 m_numDEdxBins, 0, m_dedxMaxPossible);
571 // the pion payload
572 TH2F* hDstarPiBetaGamma = new TH2F("hist_d1_211_truncBetaGamma", "hist_d1_211_truncBetaGamma;#beta*#gamma;dEdx [arb. units]",
573 m_numBGBins,
574 binsMinEdgesPi,
575 m_numDEdxBins, 0, m_dedxMaxPossible);
576
577 inputTree->Draw(Form("KaonSVDdEdx:KaonSVDdEdxTrackMomentum/%f>>hist_d1_321_truncBetaGamma", m_KaonPDGMass),
578 "nSignalDstar_sw * (KaonSVDdEdx>0)", "goff");
579 // the pion one will be built from both pions in the Dstar decay tree
580 TH2F* hDstarPiPart1BetaGamma = static_cast<TH2F*>(hDstarPiBetaGamma->Clone("hist_d1_211_truncPart1BetaGamma"));
581 TH2F* hDstarPiPart2BetaGamma = static_cast<TH2F*>(hDstarPiBetaGamma->Clone("hist_d1_211_truncPart2BetaGamma"));
582
583 inputTree->Draw(Form("PionDSVDdEdx:PionDSVDdEdxTrackMomentum/%f>>hist_d1_211_truncPart1BetaGamma", m_PionPDGMass),
584 "nSignalDstar_sw * (PionDSVDdEdx>0)",
585 "goff");
586 inputTree->Draw(Form("SlowPionSVDdEdx:SlowPionSVDdEdxTrackMomentum/%f>>hist_d1_211_truncPart2BetaGamma", m_PionPDGMass),
587 "nSignalDstar_sw * (SlowPionSVDdEdx>0)", "goff");
588 hDstarPiBetaGamma->Add(hDstarPiPart1BetaGamma);
589 hDstarPiBetaGamma->Add(hDstarPiPart2BetaGamma);
590
591
592
593 // produce the 1D profiles
594
595
596 TH1D* PionProfileMomentum = static_cast<TH1D*>(hDstarPiMomentum->ProfileX("PionProfileMomentum"));
597 PionProfileMomentum->SetTitle("PionProfile");
598 PionProfileMomentum->GetYaxis()->SetRangeUser(0, m_dedxCutoff);
599 PionProfileMomentum->GetXaxis()->SetTitle("Momentum, GeV/c");
600 PionProfileMomentum->GetYaxis()->SetTitle("dE/dx");
601 PionProfileMomentum->SetLineColor(kRed);
602
603 TH1D* PionProfileBetaGamma = static_cast<TH1D*>(hDstarPiBetaGamma->ProfileX("PionProfileBetaGamma"));
604 if (m_CustomProfile) {
605 PionProfileBetaGamma = PrepareProfile(hDstarPiBetaGamma, "PionProfileBetaGamma");
606 }
607 PionProfileBetaGamma->SetTitle("PionProfile");
608 PionProfileBetaGamma->GetYaxis()->SetRangeUser(0, m_dedxCutoff);
609 PionProfileBetaGamma->GetXaxis()->SetTitle("#beta*#gamma");
610 PionProfileBetaGamma->GetYaxis()->SetTitle("dE/dx");
611 PionProfileBetaGamma->SetLineColor(kRed);
612
613
614 TH1D* KaonProfileMomentum = static_cast<TH1D*>(hDstarKMomentum->ProfileX("KaonProfileMomentum"));
615 KaonProfileMomentum->SetTitle("KaonProfile");
616 KaonProfileMomentum->GetYaxis()->SetRangeUser(0, m_dedxCutoff);
617 KaonProfileMomentum->GetXaxis()->SetTitle("Momentum, GeV/c");
618 KaonProfileMomentum->GetYaxis()->SetTitle("dE/dx");
619 KaonProfileMomentum->SetLineColor(kRed);
620
621
622 TH1D* KaonProfileBetaGamma = static_cast<TH1D*>(hDstarKBetaGamma->ProfileX("KaonProfileBetaGamma"));
623 if (m_CustomProfile) {
624 KaonProfileBetaGamma = PrepareProfile(hDstarKBetaGamma, "KaonProfileBetaGamma");
625 }
626 KaonProfileBetaGamma->SetTitle("KaonProfile");
627
628 KaonProfileBetaGamma->GetYaxis()->SetRangeUser(0, m_dedxCutoff);
629 KaonProfileBetaGamma->GetXaxis()->SetTitle("#beta*#gamma");
630 KaonProfileBetaGamma->GetYaxis()->SetTitle("dE/dx");
631 KaonProfileBetaGamma->SetLineColor(kRed);
632
633
634
635 // normalisation
636 hDstarKMomentum->Sumw2();
637 hDstarKMomentum = Normalise2DHisto(hDstarKMomentum);
638
639 hDstarPiMomentum->Sumw2();
640 hDstarPiMomentum = Normalise2DHisto(hDstarPiMomentum);
641
642 std::unique_ptr<TList> histList(new TList);
643 histList->Add(KaonProfileMomentum);
644 histList->Add(KaonProfileBetaGamma);
645 histList->Add(hDstarKMomentum);
646
647 histList->Add(PionProfileMomentum);
648 histList->Add(PionProfileBetaGamma);
649 histList->Add(hDstarPiMomentum);
650
651 if (m_isMakePlots) {
652 TFile DstarHistogrammingPlotFile("SVDdEdxCalibrationDstarHistogramming.root", "RECREATE");
653 histList->Write();
654 DstarHistogrammingPlotFile.Close();
655 }
656
657 return histList;
658}
659
660
661std::unique_ptr<TList> SVDdEdxCalibrationAlgorithm::GammaHistogramming(std::shared_ptr<TTree> preselTree)
662{
663 B2INFO("Histogramming the converted photon selection...");
664 gROOT->SetBatch(true);
665
666
667 if (preselTree->GetEntries() == 0) {
668 B2FATAL("The Gamma tree is empty, stopping here");
669 }
670 preselTree->SetEstimate(-1);
671 std::vector<double> pbins = CreatePBinningScheme();
672
673
674 TH2F* hGammaEMomentum = new TH2F("hist_d1_11_truncMomentum", "hist_d1_11_trunc;Momentum [GeV/c];dEdx [arb. units]", m_numPBins,
675 pbins.data(), m_numDEdxBins, 0, m_dedxCutoff);
676
677 TH2F* hGammaEPart1Momentum = static_cast<TH2F*>(hGammaEMomentum->Clone("hist_d1_11_truncPart1Momentum"));
678 TH2F* hGammaEPart2Momentum = static_cast<TH2F*>(hGammaEMomentum->Clone("hist_d1_11_truncPart2Momentum"));
679
680 preselTree->Draw("FirstElectronSVDdEdx:FirstElectronSVDdEdxTrackMomentum>>hist_d1_11_truncPart1Momentum",
681 "FirstElectronSVDdEdx>0 && DIRA>0.995 && dr>1.2", "goff");
682 preselTree->Draw("SecondElectronSVDdEdx:SecondElectronSVDdEdxTrackMomentum>>hist_d1_11_truncPart2Momentum",
683 "SecondElectronSVDdEdx>0 && DIRA>0.995 && dr>1.2", "goff");
684 hGammaEMomentum->Add(hGammaEPart1Momentum);
685 hGammaEMomentum->Add(hGammaEPart2Momentum);
686
687// get a distribution that contains both pions to get a typical kinematics for the binning scheme
688 preselTree->Draw(
689 Form("FirstElectronSVDdEdxTrackMomentum/%f* (event %% 2==0) + SecondElectronSVDdEdxTrackMomentum/%f* (event %% 2==1)",
690 m_ElectronPDGMass, m_ElectronPDGMass),
691 "", "goff", ((preselTree->GetEntries()) / m_numBGBins)*m_numBGBins);
692 double* ElectronMomentumDataset = preselTree->GetV1();
693
694 TKDTreeBinning* kdBinsE = new TKDTreeBinning(((preselTree->GetEntries()) / m_numBGBins)*m_numBGBins, 1, ElectronMomentumDataset,
695 m_numBGBins);
696 const double* binsMinEdgesEOriginal = kdBinsE->SortOneDimBinEdges();
697 double* binsMinEdgesE = const_cast<double*>(binsMinEdgesEOriginal);
698 binsMinEdgesE[0] = 0.;
699 binsMinEdgesE[m_numBGBins + 1] = 10000.;
700
701
702 TH2F* hGammaEBetaGamma = new TH2F("hist_d1_11_truncBetaGamma", "hist_d1_11_truncBetaGamma;#beta*#gamma;dEdx [arb. units]",
703 m_numBGBins,
704 binsMinEdgesE,
705 m_numDEdxBins, 0, m_dedxMaxPossible);
706 TH2F* hGammaEPart1BetaGamma = static_cast<TH2F*>(hGammaEBetaGamma->Clone("hist_d1_11_truncPart1BetaGamma"));
707 TH2F* hGammaEPart2BetaGamma = static_cast<TH2F*>(hGammaEBetaGamma->Clone("hist_d1_11_truncPart2BetaGamma"));
708
709 preselTree->Draw(Form("FirstElectronSVDdEdx:FirstElectronSVDdEdxTrackMomentum/%f>>hist_d1_11_truncPart1BetaGamma",
710 m_ElectronPDGMass),
711 "FirstElectronSVDdEdx>0 && DIRA>0.995 && dr>1.2 && FirstElectronSVDdEdx<1.8e6", "goff");
712 preselTree->Draw(Form("SecondElectronSVDdEdx:SecondElectronSVDdEdxTrackMomentum/%f>>hist_d1_11_truncPart2BetaGamma",
713 m_ElectronPDGMass),
714 "SecondElectronSVDdEdx>0 && DIRA>0.995 && dr>1.2 && SecondElectronSVDdEdx<1.8e6", "goff");
715 hGammaEBetaGamma->Add(hGammaEPart1BetaGamma);
716 hGammaEBetaGamma->Add(hGammaEPart2BetaGamma);
717
718
719 // produce the 1D profile (for data-MC comparisons)
720 TH1D* ElectronProfileMomentum = static_cast<TH1D*>(hGammaEMomentum->ProfileX("ElectronProfileMomentum"));
721 ElectronProfileMomentum->SetTitle("ElectronProfile");
722 ElectronProfileMomentum->GetYaxis()->SetRangeUser(0, m_dedxCutoff);
723 ElectronProfileMomentum->GetXaxis()->SetTitle("Momentum, GeV/c");
724 ElectronProfileMomentum->GetYaxis()->SetTitle("dE/dx");
725 ElectronProfileMomentum->SetLineColor(kRed);
726
727// beta*gamma profile for the fit
728 TH1D* ElectronProfileBetaGamma = static_cast<TH1D*>(hGammaEBetaGamma->ProfileX("ElectronProfileBetaGamma"));
729 if (m_CustomProfile) {
730 ElectronProfileBetaGamma = PrepareProfile(hGammaEBetaGamma, "ElectronProfileBetaGamma");
731 }
732 ElectronProfileBetaGamma->SetTitle("ElectronProfile");
733
734 ElectronProfileBetaGamma->GetYaxis()->SetRangeUser(0, m_dedxCutoff);
735 ElectronProfileBetaGamma->GetXaxis()->SetTitle("#beta*#gamma");
736 ElectronProfileBetaGamma->GetYaxis()->SetTitle("dE/dx");
737 ElectronProfileBetaGamma->SetLineColor(kRed);
738
739
740 hGammaEMomentum = Normalise2DHisto(hGammaEMomentum);
741
742 std::unique_ptr<TList> histList(new TList);
743 histList->Add(ElectronProfileMomentum);
744 histList->Add(ElectronProfileBetaGamma);
745 histList->Add(hGammaEMomentum);
746
747 if (m_isMakePlots) {
748 TFile GammaHistogrammingPlotFile("SVDdEdxCalibrationGammaHistogramming.root", "RECREATE");
749 histList->Write();
750 GammaHistogrammingPlotFile.Close();
751 }
752
753 return histList;
754
755}
756
757std::unique_ptr<TList> SVDdEdxCalibrationAlgorithm::GenerateNewHistograms(std::shared_ptr<TTree> ttreeLambda,
758 std::shared_ptr<TTree> ttreeDstar,
759 std::shared_ptr<TTree> ttreeGamma, std::shared_ptr<TTree> ttreeGeneric)
760{
761 gROOT->SetBatch(true);
762 gStyle->SetOptStat(0);
763// run the background subtraction and histogramming parts
764 TTree* treeLambda = LambdaMassFit(ttreeLambda);
765 std::unique_ptr<TList> HistListLambda = LambdaHistogramming(treeLambda);
766 TH1D* ProtonProfileBetaGamma = static_cast<TH1D*>(HistListLambda->FindObject("ProtonProfileBetaGamma"));
767 TH2F* Proton2DHistogram = static_cast<TH2F*>(HistListLambda->FindObject("hist_d1_2212_truncMomentum"));
768
769 TTree* treeDstar = DstarMassFit(ttreeDstar);
770 std::unique_ptr<TList> HistListDstar = DstarHistogramming(treeDstar);
771 TH1D* PionProfileBetaGamma = static_cast<TH1D*>(HistListDstar->FindObject("PionProfileBetaGamma"));
772 TH2F* Pion2DHistogram = static_cast<TH2F*>(HistListDstar->FindObject("hist_d1_211_truncMomentum"));
773 TH1D* KaonProfileBetaGamma = static_cast<TH1D*>(HistListDstar->FindObject("KaonProfileBetaGamma"));
774 TH2F* Kaon2DHistogram = static_cast<TH2F*>(HistListDstar->FindObject("hist_d1_321_truncMomentum"));
775
776 std::unique_ptr<TList> HistListGamma = GammaHistogramming(ttreeGamma);
777 TH1D* ElectronProfileBetaGamma = static_cast<TH1D*>(HistListGamma->FindObject("ElectronProfileBetaGamma"));
778 TH2F* Electron2DHistogram = static_cast<TH2F*>(HistListGamma->FindObject("hist_d1_11_truncMomentum"));
779
780 int cred = TColor::GetColor("#e31a1c");
781 PionProfileBetaGamma->SetMarkerSize(4);
782 PionProfileBetaGamma->SetLineWidth(2);
783 PionProfileBetaGamma->SetMarkerColor(cred);
784 PionProfileBetaGamma->SetLineColor(cred);
785
786 int cpink = TColor::GetColor("#807dba");
787 KaonProfileBetaGamma->SetMarkerSize(4);
788 KaonProfileBetaGamma->SetLineWidth(2);
789 KaonProfileBetaGamma->SetMarkerColor(cpink);
790 KaonProfileBetaGamma->SetLineColor(cpink);
791
792 int cblue = TColor::GetColor("#084594");
793 ProtonProfileBetaGamma->SetMarkerSize(4);
794 ProtonProfileBetaGamma->SetLineWidth(2);
795 ProtonProfileBetaGamma->SetMarkerColor(cblue);
796 ProtonProfileBetaGamma->SetLineColor(cblue);
797
798 int cgreen = TColor::GetColor("#238b45");
799 ElectronProfileBetaGamma->SetMarkerSize(4);
800 ElectronProfileBetaGamma->SetLineWidth(2);
801 ElectronProfileBetaGamma->SetMarkerColor(cgreen);
802 ElectronProfileBetaGamma->SetLineColor(cgreen);
803
804// prepare the fitting
805
806 PionProfileBetaGamma->GetYaxis()->SetRangeUser(5.e5, 5.5e6);
807 KaonProfileBetaGamma->GetYaxis()->SetRangeUser(5.e5, 5.5e6);
808 ProtonProfileBetaGamma->GetYaxis()->SetRangeUser(5.e5, 5.5e6);
809
810// enhance the proton histogram (which ends at beta*gamma around 3) by adding pion data above this value – otherwise the fit is very unstable
811 auto PionEdges = PionProfileBetaGamma->GetXaxis()->GetXbins()->GetArray();
812 auto ProtonEdges = ProtonProfileBetaGamma->GetXaxis()->GetXbins()->GetArray();
813
814 std::vector<float> CombinedEdgesVector;
815
816 double borderline = 3.;
817
818 for (int i = 0; i < ProtonProfileBetaGamma->GetNbinsX() + 1; i++)
819 if (ProtonEdges[i] < borderline) CombinedEdgesVector.push_back(ProtonEdges[i]);
820
821
822 for (int i = 0; i < PionProfileBetaGamma->GetNbinsX() + 1; i++)
823 if (PionEdges[i] > borderline) CombinedEdgesVector.push_back(PionEdges[i]);
824
825
826 TH1D* CombinedHistogramPAndPi = new TH1D("CombinedHistogramPAndPi", "histo_for_fit", CombinedEdgesVector.size() - 1,
827 CombinedEdgesVector.data());
828
829 int iterator = 1;
830 for (int i = 1; i < ProtonProfileBetaGamma->GetNbinsX() + 1; i++)
831 if (ProtonEdges[i - 1] < borderline) {
832 CombinedHistogramPAndPi->SetBinContent(i, ProtonProfileBetaGamma->GetBinContent(i));
833 CombinedHistogramPAndPi->SetBinError(i, ProtonProfileBetaGamma->GetBinError(i));
834 iterator++;
835 }
836
837 for (int i = 1; i < PionProfileBetaGamma->GetNbinsX() + 1; i++)
838 if (PionEdges[i - 1] > borderline) {
839
840 CombinedHistogramPAndPi->SetBinContent(iterator, PionProfileBetaGamma->GetBinContent(i));
841 CombinedHistogramPAndPi->SetBinError(iterator, PionProfileBetaGamma->GetBinError(i));
842 iterator++;
843 }
844
845// define the beta*gamma vs momentum function
846 TF1* BetaGammaFunctionPion = new TF1("BetaGammaFunctionPion", "[0] + [1] * x/[2] + [5]/(x^2/[2]^2 + [3])**[4] + [6]* (x/[2])**0.5",
847 0.01, 25.);
848
849 BetaGammaFunctionPion->SetNpx(1000);
850
851 BetaGammaFunctionPion->SetParameters(5.e5, 2.e3, 1, 0.15, 1.2, 6.e5, 3.e5);
852
853 BetaGammaFunctionPion->SetParLimits(0, 3.e5, 7.e5);
854 BetaGammaFunctionPion->SetParLimits(1, -3.e4, 1.e4);
855 BetaGammaFunctionPion->SetParLimits(3, 0.1, 0.2);
856 BetaGammaFunctionPion->SetParLimits(4, 0.9, 1.6);
857 BetaGammaFunctionPion->SetParLimits(5, 3.e5, 7.e5);
858 BetaGammaFunctionPion->SetParLimits(6, 0., 1.e6);
859 BetaGammaFunctionPion->FixParameter(2, 1);
860 if (m_FixUnstableFitParameter) BetaGammaFunctionPion->FixParameter(3, 0.15);
861
862// fit it to the pion data
863 ROOT::Math::MinimizerOptions::SetDefaultMinimizer("Minuit2", "Migrad");
864 auto FitResultBetaGammaPion = PionProfileBetaGamma->Fit("BetaGammaFunctionPion", "0SI", "", 0.4, 25);
865
866 if ((FitResultBetaGammaPion->Status() > 1) || (BetaGammaFunctionPion->Eval(1) < 5.e5) || (BetaGammaFunctionPion->Eval(1) > 5.e6)) {
867 BetaGammaFunctionPion->FixParameter(3, 0.15);
868 FitResultBetaGammaPion = PionProfileBetaGamma->Fit("BetaGammaFunctionPion", "0SI", "", 0.4, 25);
869 }
870 if ((FitResultBetaGammaPion->Status() > 1) || (BetaGammaFunctionPion->Eval(1) < 5.e5) || (BetaGammaFunctionPion->Eval(1) > 5.e6)) {
871 BetaGammaFunctionPion->FixParameter(3, 0.15);
872 FitResultBetaGammaPion = PionProfileBetaGamma->Fit("BetaGammaFunctionPion", "0SI", "", 0.45, 25);
873 }
874 if ((FitResultBetaGammaPion->Status() > 1) || (BetaGammaFunctionPion->Eval(1) < 5.e5) || (BetaGammaFunctionPion->Eval(1) > 5.e6)) {
875 BetaGammaFunctionPion->FixParameter(3, 0.15);
876 FitResultBetaGammaPion = PionProfileBetaGamma->Fit("BetaGammaFunctionPion", "0S", "", 0.5, 25);
877 }
878
879 B2INFO("BetaGamma fit for pions done. Fit status: " << FitResultBetaGammaPion->Status());
880 // B2INFO(FitResultBetaGammaPion->Print(Belle2::LogConfig::c_Info));
881 B2INFO("Fit parameters:");
882 B2INFO("p0: " << BetaGammaFunctionPion->GetParameter(0) << " +- " << BetaGammaFunctionPion->GetParError(0));
883 B2INFO("p1: " << BetaGammaFunctionPion->GetParameter(1) << " +- " << BetaGammaFunctionPion->GetParError(1));
884 B2INFO("p2: " << BetaGammaFunctionPion->GetParameter(2) << " +- " << BetaGammaFunctionPion->GetParError(2));
885 B2INFO("p3: " << BetaGammaFunctionPion->GetParameter(3) << " +- " << BetaGammaFunctionPion->GetParError(3));
886 B2INFO("p4: " << BetaGammaFunctionPion->GetParameter(4) << " +- " << BetaGammaFunctionPion->GetParError(4));
887 B2INFO("p5: " << BetaGammaFunctionPion->GetParameter(5) << " +- " << BetaGammaFunctionPion->GetParError(5));
888 B2INFO("p6: " << BetaGammaFunctionPion->GetParameter(6) << " +- " << BetaGammaFunctionPion->GetParError(6));
889
890 // repeat the same for kaons
891 TF1* BetaGammaFunctionKaon = new TF1("BetaGammaFunctionKaon", "[0] + [1] * x/[2] + [5]/(x^2/[2]^2 + [3])**[4]+ [6]* (x/[2])**0.5",
892 0.01, 25.);
893
894 BetaGammaFunctionKaon->SetNpx(1000);
895 BetaGammaFunctionKaon->SetParameters(5.e5, 2.e3, 1, 0.15, 1.2, 6.e5, 3.e5);
896
897 BetaGammaFunctionKaon->SetParLimits(0, 3.e5, 7.e5);
898 BetaGammaFunctionKaon->SetParLimits(1, -3.e4, 1.e4);
899 BetaGammaFunctionKaon->SetParLimits(3, 0.1, 0.2);
900 BetaGammaFunctionKaon->SetParLimits(4, 0.9, 1.6);
901 BetaGammaFunctionKaon->SetParLimits(5, 3.e5, 7.e5);
902 BetaGammaFunctionKaon->SetParLimits(6, 0., 1.e6);
903
904 BetaGammaFunctionKaon->FixParameter(2, 1);
905 if (m_FixUnstableFitParameter) BetaGammaFunctionKaon->FixParameter(3, 0.15);
906
907 BetaGammaFunctionKaon->SetLineColor(KaonProfileBetaGamma->GetMarkerColor());
908
909 auto FitResultBetaGammaKaon = KaonProfileBetaGamma->Fit("BetaGammaFunctionKaon", "0SI", "", 0.4, 8.5);
910
911 if ((FitResultBetaGammaKaon->Status() > 1) || (BetaGammaFunctionKaon->Eval(1) < 5.e5) || (BetaGammaFunctionKaon->Eval(1) > 5.e6)) {
912 BetaGammaFunctionKaon->FixParameter(3, 0.15);
913 FitResultBetaGammaKaon = KaonProfileBetaGamma->Fit("BetaGammaFunctionKaon", "0SI", "", 0.4, 8.5);
914 }
915 if ((FitResultBetaGammaKaon->Status() > 1) || (BetaGammaFunctionKaon->Eval(1) < 5.e5) || (BetaGammaFunctionKaon->Eval(1) > 5.e6)) {
916 BetaGammaFunctionKaon->FixParameter(3, 0.15);
917 FitResultBetaGammaKaon = KaonProfileBetaGamma->Fit("BetaGammaFunctionKaon", "0SI", "", 0.45, 8);
918 }
919 if ((FitResultBetaGammaKaon->Status() > 1) || (BetaGammaFunctionKaon->Eval(1) < 5.e5) || (BetaGammaFunctionKaon->Eval(1) > 5.e6)) {
920 BetaGammaFunctionKaon->FixParameter(3, 0.15);
921 FitResultBetaGammaKaon = KaonProfileBetaGamma->Fit("BetaGammaFunctionKaon", "0S", "", 0.5, 8);
922 }
923
924 B2INFO("BetaGamma fit for kaons done. Fit status: " << FitResultBetaGammaKaon->Status());
925 B2INFO("Fit parameters:");
926 B2INFO("p0: " << BetaGammaFunctionKaon->GetParameter(0) << " +- " << BetaGammaFunctionKaon->GetParError(0));
927 B2INFO("p1: " << BetaGammaFunctionKaon->GetParameter(1) << " +- " << BetaGammaFunctionKaon->GetParError(1));
928 B2INFO("p2: " << BetaGammaFunctionKaon->GetParameter(2) << " +- " << BetaGammaFunctionKaon->GetParError(2));
929 B2INFO("p3: " << BetaGammaFunctionKaon->GetParameter(3) << " +- " << BetaGammaFunctionKaon->GetParError(3));
930 B2INFO("p4: " << BetaGammaFunctionKaon->GetParameter(4) << " +- " << BetaGammaFunctionKaon->GetParError(4));
931 B2INFO("p5: " << BetaGammaFunctionKaon->GetParameter(5) << " +- " << BetaGammaFunctionKaon->GetParError(5));
932 B2INFO("p6: " << BetaGammaFunctionKaon->GetParameter(6) << " +- " << BetaGammaFunctionKaon->GetParError(6));
933
934 // repeat the same for protons
935 TF1* BetaGammaFunctionProton = new TF1("BetaGammaFunctionProton",
936 "[0] + [1] * x/[2] + [5]/(x^2/[2]^2 + [3])**[4]+ [6]* (x/[2])**0.5", 0.01, 25.);
937
938 BetaGammaFunctionProton->SetNpx(1000);
939
940 BetaGammaFunctionProton->SetParameters(5.e5, 2.e3, 1, 0.15, 1.2, 6.e5, 3.e5);
941
942 BetaGammaFunctionProton->SetParLimits(0, 3.e5, 7.e5);
943 BetaGammaFunctionProton->SetParLimits(1, -3.e4, 1.e4);
944 BetaGammaFunctionProton->SetParLimits(3, 0.1, 0.2);
945 BetaGammaFunctionProton->SetParLimits(4, 0.9, 1.6);
946 BetaGammaFunctionProton->SetParLimits(5, 3.e5, 7.e5);
947 BetaGammaFunctionProton->SetParLimits(6, 0., 1.e6);
948
949 BetaGammaFunctionProton->FixParameter(2, 1);
950 if (m_FixUnstableFitParameter) BetaGammaFunctionProton->FixParameter(3, 0.15);
951
952 BetaGammaFunctionProton->SetLineColor(ProtonProfileBetaGamma->GetMarkerColor());
953
954 auto FitResultBetaGammaProton = CombinedHistogramPAndPi->Fit("BetaGammaFunctionProton", "0SI", "", 0.45, 15);
955
956 if ((FitResultBetaGammaProton->Status() > 1) || (BetaGammaFunctionProton->Eval(1) < 5.e5)
957 || (BetaGammaFunctionProton->Eval(1) > 5.e6)) {
958 BetaGammaFunctionProton->FixParameter(3, 0.15);
959 FitResultBetaGammaProton = CombinedHistogramPAndPi->Fit("BetaGammaFunctionProton", "0SI", "", 0.45, 15);
960 }
961 if ((FitResultBetaGammaProton->Status() > 1) || (BetaGammaFunctionProton->Eval(1) < 5.e5)
962 || (BetaGammaFunctionProton->Eval(1) > 5.e6)) {
963 BetaGammaFunctionProton->FixParameter(3, 0.15);
964 FitResultBetaGammaProton = CombinedHistogramPAndPi->Fit("BetaGammaFunctionProton", "0SI", "", 0.45, 10);
965 }
966 if ((FitResultBetaGammaProton->Status() > 1) || (BetaGammaFunctionProton->Eval(1) < 5.e5)
967 || (BetaGammaFunctionProton->Eval(1) > 5.e6)) {
968 BetaGammaFunctionProton->FixParameter(3, 0.15);
969 FitResultBetaGammaProton = CombinedHistogramPAndPi->Fit("BetaGammaFunctionProton", "0S", "", 0.5, 10);
970 }
971
972 B2INFO("BetaGamma fit for protons done. Fit status: " << FitResultBetaGammaProton->Status());
973 B2INFO("Fit parameters:");
974 B2INFO("p0: " << BetaGammaFunctionProton->GetParameter(0) << " +- " << BetaGammaFunctionProton->GetParError(0));
975 B2INFO("p1: " << BetaGammaFunctionProton->GetParameter(1) << " +- " << BetaGammaFunctionProton->GetParError(1));
976 B2INFO("p2: " << BetaGammaFunctionProton->GetParameter(2) << " +- " << BetaGammaFunctionProton->GetParError(2));
977 B2INFO("p3: " << BetaGammaFunctionProton->GetParameter(3) << " +- " << BetaGammaFunctionProton->GetParError(3));
978 B2INFO("p4: " << BetaGammaFunctionProton->GetParameter(4) << " +- " << BetaGammaFunctionProton->GetParError(4));
979 B2INFO("p5: " << BetaGammaFunctionProton->GetParameter(5) << " +- " << BetaGammaFunctionProton->GetParError(5));
980 B2INFO("p6: " << BetaGammaFunctionProton->GetParameter(6) << " +- " << BetaGammaFunctionProton->GetParError(6));
981
982 if (m_isMakePlots) {
983// plot a summary of all beta*gamma fits for hadrons
984 std::unique_ptr<TCanvas> CombinedCanvasHadrons(new TCanvas("CombinedCanvasHadrons", "Hadron beta*gamma fits", 10, 10, 1000, 700));
985 gStyle->SetOptFit(1111);
986
987 PionProfileBetaGamma->Draw();
988 PionProfileBetaGamma->GetListOfFunctions()->Add(BetaGammaFunctionPion);
989 KaonProfileBetaGamma->Draw("SAME");
990 KaonProfileBetaGamma->GetListOfFunctions()->Add(BetaGammaFunctionKaon);
991 ProtonProfileBetaGamma->Draw("SAME");
992 ProtonProfileBetaGamma->GetListOfFunctions()->Add(BetaGammaFunctionProton);
993// BetaGammaFunctionPion->Draw("SAME");
994// BetaGammaFunctionKaon->Draw("SAME");
995// BetaGammaFunctionProton->Draw("SAME");
996 auto legend = new TLegend(0.4, 0.7, 0.65, 0.9);
997 legend->AddEntry(PionProfileBetaGamma, "Pions", "lep");
998 legend->AddEntry(KaonProfileBetaGamma, "Kaons", "lep");
999 legend->AddEntry(ProtonProfileBetaGamma, "Protons", "lep");
1000 legend->Draw();
1001
1002 gPad->SetLogx();
1003 gPad->SetLogy();
1004
1005 CombinedCanvasHadrons->Print("HadronBetaGammaFits.pdf");
1006 TFile HadronFitPlotFile("SVDdEdxCalibrationHadronFitPlotFile.root", "RECREATE");
1007 PionProfileBetaGamma->Write();
1008 KaonProfileBetaGamma->Write();
1009 ProtonProfileBetaGamma->Write();
1010 CombinedHistogramPAndPi->Write();
1011 BetaGammaFunctionPion->Write();
1012 BetaGammaFunctionKaon->Write();
1013 BetaGammaFunctionProton->Write();
1014 CombinedCanvasHadrons->Write();
1015 HadronFitPlotFile.Close();
1016 }
1017
1018
1019 // in case we assume that all hadrons are equal
1020 if (m_UsePionBGFunctionForEverything) {
1021 BetaGammaFunctionKaon = static_cast<TF1*>(BetaGammaFunctionPion->Clone("BetaGammaFunctionKaon"));
1022 BetaGammaFunctionProton = static_cast<TF1*>(BetaGammaFunctionPion->Clone("BetaGammaFunctionProton"));
1023 }
1024 if (m_UseProtonBGFunctionForEverything) {
1025 BetaGammaFunctionKaon = static_cast<TF1*>(BetaGammaFunctionProton->Clone("BetaGammaFunctionKaon"));
1026 BetaGammaFunctionPion = static_cast<TF1*>(BetaGammaFunctionProton->Clone("BetaGammaFunctionPion"));
1027 }
1028
1029 // sanity checks: are all fits ok?
1030 if ((FitResultBetaGammaProton->Status() > 1) || (BetaGammaFunctionProton->Eval(1) < 5.e5)
1031 || (BetaGammaFunctionProton->Eval(1) > 5.e6)) {
1032 if (FitResultBetaGammaPion->Status() == 0) {
1033 BetaGammaFunctionProton = static_cast<TF1*>(BetaGammaFunctionPion->Clone("BetaGammaFunctionProton"));
1034 } else if (FitResultBetaGammaKaon->Status() == 0) {
1035 BetaGammaFunctionProton = static_cast<TF1*>(BetaGammaFunctionKaon->Clone("BetaGammaFunctionProton"));
1036 } else {
1037 B2WARNING("Problem with the beta*gamma fit for protons, reverting to the default values");
1038 BetaGammaFunctionProton->SetParameters(450258, -10900.8, 1, 0.126797, 1.155, 641907, 86304.5);
1039 }
1040 }
1041
1042 if ((FitResultBetaGammaKaon->Status() > 1) || (BetaGammaFunctionKaon->Eval(1) < 5.e5) || (BetaGammaFunctionKaon->Eval(1) > 5.e6)) {
1043 if (FitResultBetaGammaProton->Status() == 0) {
1044 BetaGammaFunctionKaon = static_cast<TF1*>(BetaGammaFunctionProton->Clone("BetaGammaFunctionKaon"));
1045 } else if (FitResultBetaGammaPion->Status() == 0) {
1046 BetaGammaFunctionKaon = static_cast<TF1*>(BetaGammaFunctionPion->Clone("BetaGammaFunctionKaon"));
1047 } else {
1048 B2WARNING("Problem with the beta*gamma fit for kaons, reverting to the default values");
1049 BetaGammaFunctionKaon->SetParameters(543386, 3013.81, 1, 0.135517, 1.19742, 619509, 15484.4);
1050 }
1051 }
1052
1053 if ((FitResultBetaGammaPion->Status() > 1) || (BetaGammaFunctionPion->Eval(1) < 5.e5) || (BetaGammaFunctionPion->Eval(1) > 5.e6)) {
1054 if (FitResultBetaGammaKaon->Status() == 0) {
1055 BetaGammaFunctionPion = static_cast<TF1*>(BetaGammaFunctionKaon->Clone("BetaGammaFunctionPion"));
1056 } else if (FitResultBetaGammaProton->Status() == 0) {
1057 BetaGammaFunctionPion = static_cast<TF1*>(BetaGammaFunctionProton->Clone("BetaGammaFunctionPion"));
1058 } else {
1059 B2WARNING("Problem with the beta*gamma fit for pions, reverting to the default values");
1060 BetaGammaFunctionPion->SetParameters(537623, -1937.62, 1, 0.15292, 1.23803, 623678, 30400.9);
1061 }
1062 }
1063
1064
1065// electrons
1066
1067 TF1* BetaGammaFunctionElectron = new TF1("BetaGammaFunctionElectron", "[0] + [1]* x", 1, 10000.);
1068 BetaGammaFunctionElectron->SetParameters(6.e5, -1);
1069 BetaGammaFunctionElectron->SetParLimits(0, 3.e5, 8.e5);
1070 BetaGammaFunctionElectron->SetParLimits(1, -1.e5, 1.e5);
1071 auto FitResultBetaGammaElectron = ElectronProfileBetaGamma->Fit("BetaGammaFunctionElectron", "0SI", "", 100, 8000);
1072
1073
1074 if ((FitResultBetaGammaElectron->Status() > 1) || (BetaGammaFunctionElectron->Eval(1) < 3.e5)
1075 || (BetaGammaFunctionElectron->Eval(1) > 5.e6)) {
1076 FitResultBetaGammaElectron = ElectronProfileBetaGamma->Fit("BetaGammaFunctionElectron", "0S", "", 100, 10000);
1077 }
1078 B2INFO("BetaGamma fit for electrons done. Fit status: " << FitResultBetaGammaElectron->Status());
1079 B2INFO("Fit parameters:");
1080 B2INFO("p0: " << BetaGammaFunctionElectron->GetParameter(0) << " +- " << BetaGammaFunctionElectron->GetParError(0));
1081 B2INFO("p1: " << BetaGammaFunctionElectron->GetParameter(1) << " +- " << BetaGammaFunctionElectron->GetParError(1));
1082
1083
1084 ElectronProfileBetaGamma->SetMarkerSize(4);
1085 ElectronProfileBetaGamma->SetLineWidth(2);
1086 ElectronProfileBetaGamma->GetYaxis()->SetRangeUser(5e5, 1e6);
1087 ElectronProfileBetaGamma->GetListOfFunctions()->Add(BetaGammaFunctionElectron);
1088
1089 if (m_isMakePlots) {
1090 std::unique_ptr<TCanvas> ElectronCanvas(new TCanvas("ElectronCanvas", "Electron histogram", 10, 10, 1000, 700));
1091 ElectronProfileBetaGamma->Draw();
1092
1093 gPad->SetLogx();
1094
1095 ElectronCanvas->Print("ElectronBetaGammaFits.pdf");
1096 TFile ElectronFitPlotFile("SVDdEdxCalibrationElectronFitPlotFile.root", "RECREATE");
1097 ElectronProfileBetaGamma->Write();
1098 BetaGammaFunctionElectron->Write();
1099 ElectronCanvas->Write();
1100 ElectronFitPlotFile.Close();
1101 }
1102
1103 TF1* MomentumFunctionElectron = static_cast<TF1*>(BetaGammaFunctionElectron->Clone("MomentumFunctionElectron"));
1104 MomentumFunctionElectron->SetParameter(2, m_ElectronPDGMass);
1105 MomentumFunctionElectron->SetRange(0.01, 5.5);
1106 MomentumFunctionElectron->SetLineColor(kRed);
1107 MomentumFunctionElectron->SetLineWidth(4);
1108
1109 TF1* MomentumFunctionPion = static_cast<TF1*>(BetaGammaFunctionPion->Clone("MomentumFunctionPion"));
1110 MomentumFunctionPion->SetParameter(2, m_PionPDGMass);
1111 MomentumFunctionPion->SetRange(0.01, 5.5);
1112 MomentumFunctionPion->SetLineColor(kRed);
1113 MomentumFunctionPion->SetLineWidth(4);
1114
1115 TF1* MomentumFunctionProton = static_cast<TF1*>(BetaGammaFunctionProton->Clone("MomentumFunctionProton"));
1116 MomentumFunctionProton->SetParameter(2, m_ProtonPDGMass);
1117 MomentumFunctionProton->SetRange(0.01, 5.5);
1118 MomentumFunctionProton->SetLineColor(kRed);
1119 MomentumFunctionProton->SetLineWidth(4);
1120
1121 TF1* MomentumFunctionKaon = static_cast<TF1*>(BetaGammaFunctionKaon->Clone("MomentumFunctionKaon"));
1122 MomentumFunctionKaon->SetParameter(2, m_KaonPDGMass);
1123 MomentumFunctionKaon->SetRange(0.01, 5.5);
1124 MomentumFunctionKaon->SetLineColor(kRed);
1125 MomentumFunctionKaon->SetLineWidth(4);
1126
1127 gStyle->SetOptFit(1111);
1128 std::unique_ptr<TCanvas> CanvasOverlays(new TCanvas("CanvasOverlays", "overlays", 1300, 1000));
1129 CanvasOverlays->Divide(2, 2);
1130 CanvasOverlays->cd(1); Electron2DHistogram->Draw(); MomentumFunctionElectron->Draw("SAME");
1131 CanvasOverlays->cd(2); Pion2DHistogram->Draw(); MomentumFunctionPion->Draw("SAME");
1132 CanvasOverlays->cd(3); Kaon2DHistogram->Draw(); MomentumFunctionKaon->Draw("SAME");
1133 CanvasOverlays->cd(4); Proton2DHistogram->Draw(); MomentumFunctionProton->Draw("SAME");
1134 CanvasOverlays->Print("SVDdEdxOverlaysFitsHistos.pdf");
1135
1136 TF1* MomentumFunctionDeuteron = static_cast<TF1*>(BetaGammaFunctionProton->Clone("MomentumFunctionDeuteron"));
1137 MomentumFunctionDeuteron->SetParameter(2, m_DeuteronPDGMass);
1138 MomentumFunctionDeuteron->SetRange(0.01, 5.5);
1139 MomentumFunctionDeuteron->SetLineColor(kRed);
1140
1141 TF1* MomentumFunctionMuon = static_cast<TF1*>(BetaGammaFunctionPion->Clone("MomentumFunctionMuon"));
1142 MomentumFunctionMuon->SetParameter(2, m_MuonPDGMass);
1143 MomentumFunctionMuon->SetRange(0.01, 5.5);
1144 MomentumFunctionMuon->SetLineColor(kRed);
1145
1146// overlay all fits in one plot
1147
1148 std::unique_ptr<TCanvas> OverlayAllTracksCanvas(new TCanvas("OverlayAllTracksCanvas", "The Ultimate Plot", 10, 10, 1000, 700));
1149
1150 TH2F* AllTracksHistogram = new TH2F("AllTracksHistogram", "AllTracksHistogram;Momentum [GeV/c];dEdx [arb. units]", 1000, 0.05, 5,
1151 1000, 2.e5, 6.e6);
1152
1153 ttreeGeneric->Draw("TrackSVDdEdx:TrackSVDdEdxTrackMomentum>>AllTracksHistogram", "TracknSVDHits>7", "goff");
1154 AllTracksHistogram->Draw("COLZ");
1155 AllTracksHistogram->GetXaxis()->SetTitle("Momentum [GeV/c]");
1156 AllTracksHistogram->GetYaxis()->SetTitle("dE/dx [arbitrary units]");
1157 MomentumFunctionElectron->Draw("SAME");
1158 MomentumFunctionMuon->Draw("SAME");
1159 MomentumFunctionPion->Draw("SAME");
1160 MomentumFunctionKaon->Draw("SAME");
1161 MomentumFunctionProton->Draw("SAME");
1162 MomentumFunctionDeuteron->Draw("SAME");
1163 OverlayAllTracksCanvas->SetLogx();
1164 OverlayAllTracksCanvas->SetLogz();
1165
1166 OverlayAllTracksCanvas->Print("SVDdEdxAllTracksWithFits.pdf");
1167 TFile OverlayAllTracksPlotFile("SVDdEdxCalibrationOverlayAllTracks.root", "RECREATE");
1168 AllTracksHistogram->Write();
1169 MomentumFunctionElectron->Write();
1170 MomentumFunctionMuon->Write();
1171 MomentumFunctionPion->Write();
1172 MomentumFunctionKaon->Write();
1173 MomentumFunctionProton->Write();
1174 MomentumFunctionDeuteron->Write();
1175 OverlayAllTracksCanvas->Write();
1176 OverlayAllTracksPlotFile.Close();
1177
1178
1179// resolution studies //
1180
1181// For resolution measurement, we need to take a ProjectionY of the data histograms in the momentum range where the dEdx is flat vs momentum. We use our educated guess of the flat range (e.g. 0.6-1 GeV for pions) and FindBin to figure out which bin numbers those momentum values correspond to.
1182 double PionRangeMin = 0.6;
1183 double PionRangeMax = 1.;
1184 double KaonRangeMin = 1.9;
1185 double KaonRangeMax = 3;
1186 double ElectronRangeMin = 1.;
1187 double ElectronRangeMax = 1.4;
1188
1189 auto PionResolutionHistogram = Pion2DHistogram->ProjectionY("PionResolutionHistogram",
1190 Pion2DHistogram->GetXaxis()->FindBin(PionRangeMin),
1191 Pion2DHistogram->GetXaxis()->FindBin(PionRangeMax));
1192 auto ElectronResolutionHistogram = Electron2DHistogram->ProjectionY("ElectronResolutionHistogram",
1193 Electron2DHistogram->GetXaxis()->FindBin(ElectronRangeMin), Electron2DHistogram->GetXaxis()->FindBin(ElectronRangeMax));
1194 auto KaonResolutionHistogram = Kaon2DHistogram->ProjectionY("KaonResolutionHistogram",
1195 Kaon2DHistogram->GetXaxis()->FindBin(KaonRangeMin),
1196 Kaon2DHistogram->GetXaxis()->FindBin(KaonRangeMax));
1197// for protons, there is not enough data in the flat range.
1198
1199
1200 TF1* PionResolutionFunction = new TF1("PionResolutionFunction",
1201 "[0]*TMath::Landau(x, [1], [1]*[2])*TMath::Gaus(x, [1], [1]*[2]*[4]) + [3]*TMath::Gaus(x, [1], [1]*[2]*[5])", 100e3, 1500e3);
1202// parameter [1] is the mean of the Landau
1203// parameter [2] is the relative resolution (w.r.t. the mean) of the Landau
1204// parameters [4]-[5] are the relative resolution of Gauss contributions w.r.t. that of Landau
1205// parameters [0] and [3] are fractions of the two components
1206 PionResolutionFunction->SetParameters(1, 6.e5, 0.1, 0.5, 2, 1);
1207 PionResolutionFunction->SetParLimits(0, 0, 500);
1208 PionResolutionFunction->SetParLimits(1, 3.e5, 8.e5);
1209 PionResolutionFunction->SetParLimits(2, 0, 1);
1210 PionResolutionFunction->SetParLimits(3, 0, 500);
1211 PionResolutionFunction->SetParLimits(4, 0, 7);
1212 PionResolutionFunction->SetParLimits(5, 1, 7);
1213 PionResolutionFunction->SetNpx(1000);
1214 auto FitResultResolutionPion = PionResolutionHistogram->Fit(PionResolutionFunction, "RSI");
1215
1216 B2INFO("relative resolution for pions: " << PionResolutionFunction->GetParameter(2));
1217 B2INFO("resolution for pions: fit status" << FitResultResolutionPion->Status());
1218
1219 TF1* KaonResolutionFunction = new TF1("KaonResolutionFunction",
1220 "[0]*TMath::Landau(x, [1], [1]*[2])*TMath::Gaus(x, [1], [1]*[2]*[4]) + [3]*TMath::Gaus(x, [1], [1]*[2]*[5])", 100e3, 1500e3);
1221
1222
1223 KaonResolutionFunction->SetParameters(1, 6.e5, 0.1, 0.5, 2, 1);
1224 KaonResolutionFunction->SetParLimits(0, 0, 500);
1225 KaonResolutionFunction->SetParLimits(1, 3.e5, 8.e5);
1226 KaonResolutionFunction->SetParLimits(2, 0, 1);
1227 KaonResolutionFunction->SetParLimits(3, 0, 500);
1228 KaonResolutionFunction->SetParLimits(4, 0, 7);
1229 KaonResolutionFunction->SetParLimits(5, 1, 7);
1230 KaonResolutionFunction->SetNpx(1000);
1231 auto FitResultResolutionKaon = KaonResolutionHistogram->Fit(KaonResolutionFunction, "RSI");
1232
1233 B2INFO("relative resolution for kaons: " << KaonResolutionFunction->GetParameter(2));
1234 B2INFO("resolution for kaons: fit status" << FitResultResolutionKaon->Status());
1235
1236 if ((FitResultResolutionKaon->Status() > 1)
1237 && (FitResultResolutionPion->Status() <= 1)) KaonResolutionFunction = static_cast<TF1*>
1238 (PionResolutionFunction->Clone("KaonResolutionFunction"));
1239
1240
1241
1242
1243 TF1* ElectronResolutionFunction = new TF1("ElectronResolutionFunction",
1244 "[0]*TMath::Landau(x, [1], [1]*[2])*TMath::Gaus(x, [1], [1]*[2]*[4]) + [3]*TMath::Gaus(x, [1], [1]*[2]*[5])", 50e3, 1500e3);
1245
1246
1247 ElectronResolutionFunction->SetParameters(1, 6.e5, 0.1, 0.5, 2, 1);
1248 ElectronResolutionFunction->SetParLimits(0, 0, 500);
1249 ElectronResolutionFunction->SetParLimits(1, 3.e5, 8.e5);
1250 ElectronResolutionFunction->SetParLimits(2, 0, 1);
1251 ElectronResolutionFunction->SetParLimits(3, 0, 500);
1252 ElectronResolutionFunction->SetParLimits(4, 0, 7);
1253 ElectronResolutionFunction->SetParLimits(5, 1, 7);
1254 ElectronResolutionFunction->SetNpx(1000);
1255 auto FitResultResolutionElectron = ElectronResolutionHistogram->Fit(ElectronResolutionFunction, "RSI");
1256
1257 B2INFO("relative resolution for electrons: " << ElectronResolutionFunction->GetParameter(2));
1258 B2INFO("resolution for electrons: fit status" << FitResultResolutionElectron->Status());
1259
1260 // plot all the resolution fits
1261 if (m_isMakePlots) {
1262 TCanvas* CanvasResolutions = new TCanvas("CanvasResolutions", "Resolutions", 1200, 650);
1263 CanvasResolutions->Divide(3, 1);
1264 CanvasResolutions->cd(1); PionResolutionHistogram->Draw();
1265 CanvasResolutions->cd(2); KaonResolutionHistogram->Draw();
1266 CanvasResolutions->cd(3); ElectronResolutionHistogram->Draw();
1267
1268 CanvasResolutions->Print("SVDdEdxResolutions.pdf");
1269 TFile OverlayResolutionsPlotFile("SVDdEdxCalibrationResolutions.root", "RECREATE");
1270 PionResolutionHistogram->Write();
1271 KaonResolutionHistogram->Write();
1272 ElectronResolutionHistogram->Write();
1273 CanvasResolutions->Write();
1274 OverlayResolutionsPlotFile.Close();
1275 }
1276
1277// evaluate the bias correction:
1278// difference between the MomentumFunction prediction and the mean of the resolution function in the flat part
1279// it should be of the order -1e4, i.e. about -1% of the absolute dEdx value
1280 double BiasCorrectionPion = PionResolutionFunction->GetParameter(1) - MomentumFunctionPion->Eval((
1281 PionRangeMax + PionRangeMin) / 2.);
1282 B2INFO("BiasCorrectionPion = " << BiasCorrectionPion);
1283
1284// generate a new pion payload using the MomentumFunctionPion, PionResolutionFunction and the bias correction
1285 TH2F* Pion2DHistogramNew = PrepareNewHistogram(Pion2DHistogram, Form("%sNew", Pion2DHistogram->GetName()), MomentumFunctionPion,
1286 PionResolutionFunction, BiasCorrectionPion);
1287
1288// sanity check: residual between the generated distribution and the data one
1289 TH2F* Pion2DHistogramResidual = static_cast<TH2F*>(Pion2DHistogram->Clone("Pion2DHistogramResidual"));
1290 Pion2DHistogramResidual->Add(Pion2DHistogramNew, Pion2DHistogram, 1, -1);
1291 Pion2DHistogramResidual->SetMinimum(-0.15);
1292 Pion2DHistogramResidual->SetMaximum(0.15);
1293
1294 // repeat, for kaons
1295 double BiasCorrectionKaon = KaonResolutionFunction->GetParameter(1) - MomentumFunctionKaon->Eval((
1296 KaonRangeMax + KaonRangeMin) / 2.);
1297 B2INFO("BiasCorrectionKaon = " << BiasCorrectionKaon);
1298
1299 // for protons, we compare the flat part of the MomentumFunction (~3 GeV) with the mean of the kaon resolution function
1300 // as there's not enough stats in the flat part to extract proton resolution from data
1301 double BiasCorrectionProton = KaonResolutionFunction->GetParameter(1) - MomentumFunctionProton->Eval(3.);
1302 B2INFO("BiasCorrectionProton = " << BiasCorrectionProton);
1303
1304 if ((BiasCorrectionProton / BiasCorrectionKaon) > 1.5) BiasCorrectionProton =
1305 BiasCorrectionKaon; // probably something went wrong due to low statistics
1306
1307 // back to kaons: generate a new payload
1308 TH2F* Kaon2DHistogramNew = PrepareNewHistogram(Kaon2DHistogram, Form("%sNew", Kaon2DHistogram->GetName()), MomentumFunctionKaon,
1309 KaonResolutionFunction, BiasCorrectionKaon);
1310// residual generated - data for kaons
1311 TH2F* Kaon2DHistogramResidual = static_cast<TH2F*>(Kaon2DHistogram->Clone("Kaon2DHistogramResidual"));
1312 Kaon2DHistogramResidual->Add(Kaon2DHistogramNew, Kaon2DHistogram, 1, -1);
1313 Kaon2DHistogramResidual->SetMinimum(-0.15);
1314 Kaon2DHistogramResidual->SetMaximum(0.15);
1315
1316 // same for protons (we use the kaon resolution function as explained above)
1317 TH2F* Proton2DHistogramNew = PrepareNewHistogram(Proton2DHistogram, Form("%sNew", Proton2DHistogram->GetName()),
1318 MomentumFunctionProton,
1319 KaonResolutionFunction, BiasCorrectionProton);
1320
1321// residual for protons
1322 TH2F* Proton2DHistogramResidual = static_cast<TH2F*>(Proton2DHistogram->Clone("Proton2DHistogramResidual"));
1323 Proton2DHistogramResidual->Add(Proton2DHistogramNew, Proton2DHistogram, 1, -1);
1324 Proton2DHistogramResidual->SetMinimum(-0.15);
1325 Proton2DHistogramResidual->SetMaximum(0.15);
1326
1327// deuterons: same as protons, but use the MomentumFunctionDeuteron
1328 TH2F* Deuteron2DHistogramNew = PrepareNewHistogram(Proton2DHistogram, "Deuteron2DHistogramNew", MomentumFunctionDeuteron,
1329 KaonResolutionFunction,
1330 BiasCorrectionKaon);
1331 Deuteron2DHistogramNew->SetTitle("hist_d1_1000010020_trunc");
1332
1333// muons: same as pions, but use the MomentumFunctionMuon
1334 TH2F* Muon2DHistogramNew = PrepareNewHistogram(Pion2DHistogram, "Muon2DHistogramNew", MomentumFunctionMuon, PionResolutionFunction,
1335 BiasCorrectionPion);
1336 Muon2DHistogramNew->SetTitle("hist_d1_13_trunc");
1337
1338// same for electrons
1339 double BiasCorrectionElectron = ElectronResolutionFunction->GetParameter(1) - MomentumFunctionElectron->Eval((
1340 ElectronRangeMax + ElectronRangeMin) / 2.);
1341 B2INFO("BiasCorrectionElectron = " << BiasCorrectionElectron);
1342 TH2F* Electron2DHistogramNew = PrepareNewHistogram(Electron2DHistogram, Form("%sNew", Electron2DHistogram->GetName()),
1343 MomentumFunctionElectron,
1344 ElectronResolutionFunction, BiasCorrectionElectron);
1345
1346 TH2F* Electron2DHistogramResidual = static_cast<TH2F*>(Electron2DHistogram->Clone("Electron2DHistogramResidual"));
1347 Electron2DHistogramResidual->Add(Electron2DHistogramNew, Electron2DHistogram, 1, -1);
1348 Electron2DHistogramResidual->SetMinimum(-0.15);
1349 Electron2DHistogramResidual->SetMaximum(0.15);
1350
1351 Electron2DHistogramNew->SetName("Electron2DHistogramNew");
1352 Muon2DHistogramNew->SetName("Muon2DHistogramNew");
1353 Pion2DHistogramNew->SetName("Pion2DHistogramNew");
1354 Kaon2DHistogramNew->SetName("Kaon2DHistogramNew");
1355 Proton2DHistogramNew->SetName("Proton2DHistogramNew");
1356 Deuteron2DHistogramNew->SetName("Deuteron2DHistogramNew");
1357
1358// plot the summary of all the distributions
1359 if (m_isMakePlots) {
1360 TCanvas* CanvasSummaryGenerated = new TCanvas("CanvasSummaryGenerated", "Generated payloads", 1700, 850);
1361 CanvasSummaryGenerated->Divide(3, 2);
1362 CanvasSummaryGenerated->cd(1); Electron2DHistogramNew->Draw("COLZ");
1363 CanvasSummaryGenerated->cd(2); Muon2DHistogramNew->Draw("COLZ");
1364 CanvasSummaryGenerated->cd(3); Pion2DHistogramNew->Draw("COLZ");
1365 CanvasSummaryGenerated->cd(4); Kaon2DHistogramNew->Draw("COLZ");
1366 CanvasSummaryGenerated->cd(5); Proton2DHistogramNew->Draw("COLZ");
1367 CanvasSummaryGenerated->cd(6); Deuteron2DHistogramNew->Draw("COLZ");
1368
1369 CanvasSummaryGenerated->Print("SVDdEdxGeneratedPayloads.pdf");
1370 TFile SummaryGeneratedPlotFile("SVDdEdxCalibrationSummaryGenerated.root", "RECREATE");
1371 Electron2DHistogramNew->Write();
1372 Muon2DHistogramNew->Write();
1373 Pion2DHistogramNew->Write();
1374 Kaon2DHistogramNew->Write();
1375 Proton2DHistogramNew->Write();
1376 Deuteron2DHistogramNew->Write();
1377 SummaryGeneratedPlotFile.Close();
1378
1379
1380 TCanvas* CanvasSummaryData = new TCanvas("CanvasSummaryData", "Data distributions", 1700, 850);
1381 CanvasSummaryData->Divide(3, 2);
1382 CanvasSummaryData->cd(1); Electron2DHistogram->Draw("COLZ");
1383 CanvasSummaryData->cd(3); Pion2DHistogram->Draw("COLZ");
1384 CanvasSummaryData->cd(4); Kaon2DHistogram->Draw("COLZ");
1385 CanvasSummaryData->cd(5); Proton2DHistogram->Draw("COLZ");
1386
1387 CanvasSummaryData->Print("SVDdEdxDataDistributions.pdf");
1388 TFile SummaryDataPlotFile("SVDdEdxCalibrationSummaryData.root", "RECREATE");
1389 Electron2DHistogram->Write();
1390 Pion2DHistogram->Write();
1391 Kaon2DHistogram->Write();
1392 Proton2DHistogram->Write();
1393 SummaryDataPlotFile.Close();
1394
1395
1396 TCanvas* CanvasSummaryResiduals = new TCanvas("CanvasSummaryResiduals", "Residuals", 1700, 850);
1397 CanvasSummaryResiduals->Divide(3, 2);
1398 CanvasSummaryResiduals->cd(1); Electron2DHistogramResidual->Draw("COLZ");
1399 CanvasSummaryResiduals->cd(3); Pion2DHistogramResidual->Draw("COLZ");
1400 CanvasSummaryResiduals->cd(4); Kaon2DHistogramResidual->Draw("COLZ");
1401 CanvasSummaryResiduals->cd(5); Proton2DHistogramResidual->Draw("COLZ");
1402
1403
1404 CanvasSummaryResiduals->Print("SVDdEdxResiduals.pdf");
1405 TFile SummaryResidualsPlotFile("SVDdEdxCalibrationSummaryResiduals.root", "RECREATE");
1406 Electron2DHistogramResidual->Write();
1407 Pion2DHistogramResidual->Write();
1408 Kaon2DHistogramResidual->Write();
1409 Proton2DHistogramResidual->Write();
1410 SummaryResidualsPlotFile.Close();
1411 }
1412
1413
1414 // return all the generated payloads
1415 std::unique_ptr<TList> histList(new TList);
1416 histList->Add(Electron2DHistogramNew);
1417 histList->Add(Muon2DHistogramNew);
1418 histList->Add(Pion2DHistogramNew);
1419 histList->Add(Kaon2DHistogramNew);
1420 histList->Add(Proton2DHistogramNew);
1421 histList->Add(Deuteron2DHistogramNew);
1422
1423 return histList;
1424
1425}
Belle2::CalibrationAlgorithm::saveCalibration
void saveCalibration(TClonesArray *data, const std::string &name)
Store DBArray payload with given name with default IOV.
Definition CalibrationAlgorithm.cc:297
Belle2::CalibrationAlgorithm::setDescription
void setDescription(const std::string &description)
Set algorithm description (in constructor)
Definition CalibrationAlgorithm.h:321
Belle2::CalibrationAlgorithm::getRunList
const std::vector< Calibration::ExpRun > & getRunList() const
Get the list of runs for which calibration is called.
Definition CalibrationAlgorithm.h:266
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::CalibrationAlgorithm::CalibrationAlgorithm
CalibrationAlgorithm(const std::string &collectorModuleName)
Constructor - sets the prefix for collected objects (won't be accesses until execute(....
Definition CalibrationAlgorithm.h:140
Belle2::SVDdEdxCalibrationAlgorithm::m_numPBins
int m_numPBins
the number of momentum bins for the payloads
Definition SVDdEdxCalibrationAlgorithm.h:114
Belle2::SVDdEdxCalibrationAlgorithm::m_MuonPDGMass
const double m_MuonPDGMass
PDG mass for the muon.
Definition SVDdEdxCalibrationAlgorithm.h:131
Belle2::SVDdEdxCalibrationAlgorithm::m_ElectronPDGMass
const double m_ElectronPDGMass
PDG mass for the electron.
Definition SVDdEdxCalibrationAlgorithm.h:130
Belle2::SVDdEdxCalibrationAlgorithm::m_FixUnstableFitParameter
bool m_FixUnstableFitParameter
In the dEdx:betagamma fit, there is one free parameter that makes fit convergence poor.
Definition SVDdEdxCalibrationAlgorithm.h:127
Belle2::SVDdEdxCalibrationAlgorithm::m_DeuteronPDGMass
const double m_DeuteronPDGMass
PDG mass for the deuteron.
Definition SVDdEdxCalibrationAlgorithm.h:135
Belle2::SVDdEdxCalibrationAlgorithm::CreatePBinningScheme
std::vector< double > CreatePBinningScheme()
build the binning scheme for the momentum
Definition SVDdEdxCalibrationAlgorithm.h:140
Belle2::SVDdEdxCalibrationAlgorithm::Normalise2DHisto
TH2F * Normalise2DHisto(TH2F *HistoToNormalise)
Normalise a given dEdx:momentum histogram in each momentum bin, so that sum of entries in each moment...
Definition SVDdEdxCalibrationAlgorithm.h:163
Belle2::SVDdEdxCalibrationAlgorithm::m_UseProtonBGFunctionForEverything
bool m_UseProtonBGFunctionForEverything
Assume that the dEdx:betagamma trend is the same for all hadrons; use the proton trend as representat...
Definition SVDdEdxCalibrationAlgorithm.h:125
Belle2::SVDdEdxCalibrationAlgorithm::m_numBGBins
int m_numBGBins
the number of beta*gamma bins for the profile and fitting
Definition SVDdEdxCalibrationAlgorithm.h:115
Belle2::SVDdEdxCalibrationAlgorithm::DstarHistogramming
std::unique_ptr< TList > DstarHistogramming(TTree *inputTree)
produce histograms for K/pi
Definition SVDdEdxCalibrationAlgorithm.cc:517
Belle2::SVDdEdxCalibrationAlgorithm::m_dedxMaxPossible
double m_dedxMaxPossible
the approximate max possible value of dEdx
Definition SVDdEdxCalibrationAlgorithm.h:117
Belle2::SVDdEdxCalibrationAlgorithm::m_isMakePlots
bool m_isMakePlots
produce plots for monitoring
Definition SVDdEdxCalibrationAlgorithm.h:104
Belle2::SVDdEdxCalibrationAlgorithm::m_MinEvtsPerTree
int m_MinEvtsPerTree
number of events in TTree below which we don't try to fit
Definition SVDdEdxCalibrationAlgorithm.h:118
Belle2::SVDdEdxCalibrationAlgorithm::PrepareProfile
TH1D * PrepareProfile(TH2F *DataHistogram, TString NewName)
Reimplement the Profile histogram calculation for a 2D histogram.
Definition SVDdEdxCalibrationAlgorithm.h:232
Belle2::SVDdEdxCalibrationAlgorithm::m_numDEdxBins
int m_numDEdxBins
the number of dEdx bins for the payloads
Definition SVDdEdxCalibrationAlgorithm.h:113
Belle2::SVDdEdxCalibrationAlgorithm::LambdaMassFit
TTree * LambdaMassFit(std::shared_ptr< TTree > preselTree)
Mass fit for Lambda->ppi.
Definition SVDdEdxCalibrationAlgorithm.cc:164
Belle2::SVDdEdxCalibrationAlgorithm::LambdaHistogramming
std::unique_ptr< TList > LambdaHistogramming(TTree *inputTree)
produce histograms for protons
Definition SVDdEdxCalibrationAlgorithm.cc:308
Belle2::SVDdEdxCalibrationAlgorithm::m_KaonPDGMass
const double m_KaonPDGMass
PDG mass for the charged kaon.
Definition SVDdEdxCalibrationAlgorithm.h:133
Belle2::SVDdEdxCalibrationAlgorithm::PrepareNewHistogram
TH2F * PrepareNewHistogram(TH2F *DataHistogram, TString NewName, TF1 *betagamma_function, TF1 *ResolutionFunctionOriginal, double bias_correction)
Generate a new dEdx:momentum histogram from a function that encodes dEdx:momentum trend and a functio...
Definition SVDdEdxCalibrationAlgorithm.h:191
Belle2::SVDdEdxCalibrationAlgorithm::GammaHistogramming
std::unique_ptr< TList > GammaHistogramming(std::shared_ptr< TTree > preselTree)
produce histograms for e
Definition SVDdEdxCalibrationAlgorithm.cc:661
Belle2::SVDdEdxCalibrationAlgorithm::DstarMassFit
TTree * DstarMassFit(std::shared_ptr< TTree > preselTree)
Mass fit for D*->Dpi.
Definition SVDdEdxCalibrationAlgorithm.cc:384
Belle2::SVDdEdxCalibrationAlgorithm::m_UsePionBGFunctionForEverything
bool m_UsePionBGFunctionForEverything
Assume that the dEdx:betagamma trend is the same for all hadrons; use the pion trend as representativ...
Definition SVDdEdxCalibrationAlgorithm.h:123
Belle2::SVDdEdxCalibrationAlgorithm::GenerateNewHistograms
std::unique_ptr< TList > GenerateNewHistograms(std::shared_ptr< TTree > ttreeLambda, std::shared_ptr< TTree > ttreeDstar, std::shared_ptr< TTree > ttreeGamma, std::shared_ptr< TTree > ttreeGeneric)
generate high-statistics histograms
Definition SVDdEdxCalibrationAlgorithm.cc:757
Belle2::SVDdEdxCalibrationAlgorithm::calibrate
virtual EResult calibrate() override
run algorithm on data
Definition SVDdEdxCalibrationAlgorithm.cc:48
Belle2::SVDdEdxCalibrationAlgorithm::m_CustomProfile
bool m_CustomProfile
reimplement profile histogram calculation instead of the ROOT implementation?
Definition SVDdEdxCalibrationAlgorithm.h:122
Belle2::SVDdEdxCalibrationAlgorithm::m_dedxCutoff
double m_dedxCutoff
the upper edge of the dEdx binning for the payloads
Definition SVDdEdxCalibrationAlgorithm.h:116
Belle2::SVDdEdxCalibrationAlgorithm::m_ProtonPDGMass
const double m_ProtonPDGMass
PDG mass for the proton.
Definition SVDdEdxCalibrationAlgorithm.h:134
Belle2::SVDdEdxCalibrationAlgorithm::m_PionPDGMass
const double m_PionPDGMass
PDG mass for the charged pion.
Definition SVDdEdxCalibrationAlgorithm.h:132
Belle2::SVDdEdxPDFs
Specialized class for holding the SVD dE/dx PDFs.
Definition SVDdEdxPDFs.h:26
Belle2::CalibrationAlgorithm::getObjectPtr
std::shared_ptr< T > getObjectPtr(const std::string &name, const std::vector< Calibration::ExpRun > &requestedRuns)
Get calibration data object by name and list of runs, the Merge function will be called to generate t...
Definition CalibrationAlgorithm.h:414
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