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 = (TH2F*) GeneratedList->FindObject("Electron2DHistogramNew");
73 TH2F* histoMu = (TH2F*) GeneratedList->FindObject("Muon2DHistogramNew");
74 TH2F* histoPi = (TH2F*) GeneratedList->FindObject("Pion2DHistogramNew");
75 TH2F* histoK = (TH2F*) GeneratedList->FindObject("Kaon2DHistogramNew");
76 TH2F* histoP = (TH2F*) GeneratedList->FindObject("Proton2DHistogramNew");
77 TH2F* histoDeut = (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 RooDataSet* LambdaDatasetSWeighted = new RooDataSet(LambdaDataset->GetName(), LambdaDataset->GetTitle(), LambdaDataset,
302 *LambdaDataset->get());
303
304 TTree* treeLambdaSWeighted = LambdaDatasetSWeighted->GetClonedTree();
305 treeLambdaSWeighted->SetName("treeLambdaSWeighted");
306
307 B2INFO("Lambda: sPlot done. Proceed to histogramming");
308 return treeLambdaSWeighted;
309}
310
311std::unique_ptr<TList> SVDdEdxCalibrationAlgorithm::LambdaHistogramming(TTree* inputTree)
312{
313 gROOT->SetBatch(true);
314 inputTree->SetEstimate(-1);
315 std::vector<double> pbins = CreatePBinningScheme();
316
317 TH2F* hLambdaPMomentum = new TH2F("hist_d1_2212_truncMomentum", "hist_d1_2212_trunc;Momentum [GeV/c];dEdx [arb. units]",
318 m_numPBins, pbins.data(), m_numDEdxBins, 0,
319 m_dedxCutoff);
320
321 inputTree->Draw("ProtonSVDdEdx:ProtonSVDdEdxTrackMomentum>>hist_d1_2212_truncMomentum",
322 "nSignalLambda_sw * (ProtonSVDdEdx>0) * (ProtonSVDdEdxTrackMomentum>0.13)", "goff");
323
324// create isopopulated beta*gamma binning
325 inputTree->Draw(Form("ProtonSVDdEdxTrackMomentum/%f", m_ProtonPDGMass), "", "goff",
326 ((inputTree->GetEntries()) / m_numBGBins)*m_numBGBins);
327 double* ProtonMomentumDataset = inputTree->GetV1();
328
329 TKDTreeBinning* kdBinsP = new TKDTreeBinning(((inputTree->GetEntries()) / m_numBGBins)*m_numBGBins, 1, ProtonMomentumDataset,
330 m_numBGBins);
331 const double* binsMinEdgesP_pointer = kdBinsP->SortOneDimBinEdges();
332 double* binsMinEdgesP = const_cast<double*>(binsMinEdgesP_pointer);
333
334
335 binsMinEdgesP[0] = 0.1;
336 binsMinEdgesP[m_numBGBins + 1] = 50.;
337
338
339 TH2F* hLambdaPBetaGamma = new TH2F("hist_d1_2212_truncBetaGamma", "hist_d1_2212_truncBetaGamma;#beta*#gamma;dEdx [arb. units]",
340 m_numBGBins, binsMinEdgesP, m_numDEdxBins,
341 0,
342 m_dedxMaxPossible);
343
344 inputTree->Draw(Form("ProtonSVDdEdx:ProtonSVDdEdxTrackMomentum/%f>>hist_d1_2212_truncBetaGamma", m_ProtonPDGMass),
345 "nSignalLambda_sw * (ProtonSVDdEdx>0) * (ProtonSVDdEdxTrackMomentum>0.13) * (ProtonSVDdEdx>1.2e6 - 1.e6*ProtonSVDdEdxTrackMomentum)",
346 "goff");
347
348 // produce the 1D profile
349 // momentum: for data-MC comparisons
350
351 TH1F* ProtonProfileMomentum = (TH1F*)hLambdaPMomentum->ProfileX("ProtonProfileMomentum");
352 ProtonProfileMomentum->SetTitle("ProtonProfile");
353 ProtonProfileMomentum->GetYaxis()->SetRangeUser(0, m_dedxCutoff);
354 ProtonProfileMomentum->GetXaxis()->SetTitle("Momentum, GeV/c");
355 ProtonProfileMomentum->GetYaxis()->SetTitle("dE/dx");
356 ProtonProfileMomentum->SetLineColor(kRed);
357
358 //beta*gamma: for the fit
359 TH1F* ProtonProfileBetaGamma = (TH1F*)hLambdaPBetaGamma->ProfileX("ProtonProfileBetaGamma");
360 if (m_CustomProfile) {
361 ProtonProfileBetaGamma = PrepareProfile(hLambdaPBetaGamma, "ProtonProfileBetaGamma");
362 }
363 ProtonProfileBetaGamma->SetTitle("ProtonProfile");
364 ProtonProfileBetaGamma->GetYaxis()->SetRangeUser(0, m_dedxCutoff);
365 ProtonProfileBetaGamma->GetXaxis()->SetTitle("#beta*#gamma");
366 ProtonProfileBetaGamma->GetYaxis()->SetTitle("dE/dx");
367 ProtonProfileBetaGamma->SetLineColor(kRed);
368
369
370 // for each momentum bin, normalize the pdf
371 hLambdaPMomentum = Normalise2DHisto(hLambdaPMomentum);
372
373 std::unique_ptr<TList> histList(new TList);
374 histList->Add(ProtonProfileMomentum);
375 histList->Add(ProtonProfileBetaGamma);
376 histList->Add(hLambdaPMomentum);
377
378 if (m_isMakePlots) {
379 TFile LambdaHistogrammingPlotFile("SVDdEdxCalibrationLambdaHistogramming.root", "RECREATE");
380 histList->Write();
381 LambdaHistogrammingPlotFile.Close();
382 }
383
384 return histList;
385}
386
387TTree* SVDdEdxCalibrationAlgorithm::DstarMassFit(std::shared_ptr<TTree> preselTree)
388{
389 B2INFO("Configuring the Dstar fit...");
390 gROOT->SetBatch(true);
391 RooMsgService::instance().setGlobalKillBelow(RooFit::WARNING);
392
393 RooRealVar deltaM("deltaM", "m(D*)-m(D^{0})", 0.139545, 0.151, "GeV/c^{2}");
394
395 RooRealVar KaonMomentum("KaonMomentum", "momentum for Kaon (GeV)", -1.e8, 1.e8);
396 RooRealVar KaonSVDdEdxTrackMomentum("KaonSVDdEdxTrackMomentum", "momentum for Kaon (GeV), from the track", -1.e8, 1.e8);
397 RooRealVar KaonSVDdEdx("KaonSVDdEdx", "", -1.e8, 1.e8);
398 RooRealVar PionDMomentum("PionDMomentum", "momentum for pion (GeV)", -1.e8, 1.e8);
399 RooRealVar PionDSVDdEdxTrackMomentum("PionDSVDdEdxTrackMomentum", "momentum for pion (GeV), from the track", -1.e8, 1.e8);
400 RooRealVar PionDSVDdEdx("PionDSVDdEdx", "", -1.e8, 1.e8);
401 RooRealVar SlowPionMomentum("SlowPionMomentum", "momentum for slow pion (GeV)", -1.e8, 1.e8);
402 RooRealVar SlowPionSVDdEdxTrackMomentum("SlowPionSVDdEdxTrackMomentum", "momentum for slow pion (GeV), from the track", -1.e8,
403 1.e8);
404 RooRealVar SlowPionSVDdEdx("SlowPionSVDdEdx", "", -1.e8, 1.e8);
405
406 RooRealVar exp("exp", "experiment number", 0, 1.e5);
407 RooRealVar run("run", "run number", 0, 1.e8);
408 RooRealVar event("event", "event number", 0, 1.e10);
409
410 auto variables = new RooArgSet();
411 variables->add(deltaM);
412 variables->add(KaonMomentum);
413 variables->add(KaonSVDdEdxTrackMomentum);
414 variables->add(KaonSVDdEdx);
415 variables->add(PionDMomentum);
416 variables->add(PionDSVDdEdxTrackMomentum);
417 variables->add(PionDSVDdEdx);
418 variables->add(SlowPionMomentum);
419 variables->add(SlowPionSVDdEdxTrackMomentum);
420 variables->add(SlowPionSVDdEdx);
421 variables->add(exp);
422 variables->add(run);
423 variables->add(event);
424
425 RooDataSet* DstarDataset = new RooDataSet("DstarDataset", "DstarDataset", *variables, Import(*preselTree));
426
427 if (DstarDataset->sumEntries() == 0) {
428 B2FATAL("The Dstar dataset is empty, stopping here");
429 }
430
431 RooPlot* DstarFitFrame = DstarDataset->plotOn(deltaM.frame());
432
433 RooRealVar GaussMean("GaussMean", "GaussMean", 0.145, 0.140, 0.150);
434 RooRealVar GaussSigma1("GaussSigma1", "GaussSigma1", 0.01, 1.e-4, 1.0);
435 RooGaussian DstarGauss1("DstarGauss1", "DstarGauss1", deltaM, GaussMean, GaussSigma1);
436 RooRealVar GaussSigma2("GaussSigma2", "GaussSigma2", 0.001, 1.e-4, 1.0);
437 RooGaussian DstarGauss2("DstarGauss2", "DstarGauss2", deltaM, GaussMean, GaussSigma2);
438 RooRealVar fracGaussYield("fracGaussYield", "Fraction of two Gaussians", 0.75, 0.0, 1.0);
439 RooAddPdf DstarSignalPDF("DstarSignalPDF", "DstarGauss1+DstarGauss2", RooArgList(DstarGauss1, DstarGauss2), fracGaussYield);
440
441 RooRealVar dm0Bkg("dm0Bkg", "dm0", 0.13957018, 0.130, 0.140);
442 RooRealVar aBkg("aBkg", "a", -0.0784, -0.08, 3.0);
443 RooRealVar bBkg("bBkg", "b", -0.444713, -0.5, 0.4);
444 RooRealVar cBkg("cBkg", "c", 0.3);
445 RooDstD0BG DstarBkgPDF("DstarBkgPDF", "DstarBkgPDF", deltaM, dm0Bkg, cBkg, aBkg, bBkg);
446 RooRealVar nSignalDstar("nSignalDstar", "signal yield", 0.5 * preselTree->GetEntries(), 0, preselTree->GetEntries());
447 RooRealVar nBkgDstar("nBkgDstar", "background yield", 0.5 * preselTree->GetEntries(), 0, preselTree->GetEntries());
448 RooAddPdf totalPDFDstar("totalPDFDstar", "totalPDFDstar pdf", RooArgList(DstarSignalPDF, DstarBkgPDF),
449 RooArgList(nSignalDstar, nBkgDstar));
450
451 B2INFO("Dstar: Start fitting...");
452 RooFitResult* DstarFitResult = totalPDFDstar.fitTo(*DstarDataset, Save(kTRUE), PrintLevel(-1));
453
454 int status = DstarFitResult->status();
455 int covqual = DstarFitResult->covQual();
456 double diff = nSignalDstar.getValV() + nBkgDstar.getValV() - DstarDataset->sumEntries();
457
458 B2INFO("Dstar: Fit status: " << status << "; covariance quality: " << covqual);
459 // if the fit is not healthy, try again once before giving up, with a slightly different setup:
460 if ((status > 0) || (TMath::Abs(diff) > 1.) || (nSignalDstar.getError() < sqrt(nSignalDstar.getValV()))
461 || (nSignalDstar.getError() > (nSignalDstar.getValV()))) {
462
463 DstarFitResult = totalPDFDstar.fitTo(*DstarDataset, Save(), Strategy(2), Offset(1));
464 status = DstarFitResult->status();
465 covqual = DstarFitResult->covQual();
466 diff = nSignalDstar.getValV() + nBkgDstar.getValV() - DstarDataset->sumEntries();
467 B2INFO("Dstar: Updated fit status: " << status << "; covariance quality: " << covqual);
468 }
469
470 if ((status > 0) || (TMath::Abs(diff) > 1.) || (nSignalDstar.getError() < sqrt(nSignalDstar.getValV()))
471 || (nSignalDstar.getError() > (nSignalDstar.getValV()))) {
472 B2WARNING("Dstar: Fit problem: fit status " << status << "; sum of component yields minus the dataset yield is " << diff <<
473 "; signal yield is " << nSignalDstar.getValV() << ", while its uncertainty is " << nSignalDstar.getError());
474 }
475 if (covqual < 2) {
476 B2INFO("Dstar: Fit warning: covariance quality " << covqual);
477 }
478
479 totalPDFDstar.plotOn(DstarFitFrame, LineColor(TColor::GetColor("#4575b4")));
480
481 double chisquare = DstarFitFrame->chiSquare();
482 B2INFO("Dstar: Fit chi2 = " << chisquare);
483 totalPDFDstar.paramOn(DstarFitFrame, Layout(0.63, 0.96, 0.93), Format("NEU", AutoPrecision(2)));
484 DstarFitFrame->getAttText()->SetTextSize(0.03);
485
486 totalPDFDstar.plotOn(DstarFitFrame, Components("DstarSignalPDF"), LineColor(TColor::GetColor("#d73027")));
487 totalPDFDstar.plotOn(DstarFitFrame, Components("DstarBkgPDF"), LineColor(TColor::GetColor("#fc8d59")));
488 totalPDFDstar.plotOn(DstarFitFrame, LineColor(TColor::GetColor("#4575b4")));
489
490 DstarFitFrame->GetXaxis()->SetTitle("#Deltam [GeV/c^{2}]");
491 if (m_isMakePlots) {
492 std::unique_ptr<TCanvas> canvDstar(new TCanvas("canvDstar", "canvDstar"));
493 canvDstar->cd();
494
495 DstarFitFrame->Draw();
496
497 canvDstar->Print("SVDdEdxCalibrationFitDstar.pdf");
498 TFile DstarFitPlotFile("SVDdEdxCalibrationDstarFitPlotFile.root", "RECREATE");
499 canvDstar->Write();
500 DstarFitPlotFile.Close();
501 }
502
504
505 RooStats::SPlot* sPlotDatasetDstar = new RooStats::SPlot("sData", "An SPlot", *DstarDataset, &totalPDFDstar,
506 RooArgList(nSignalDstar, nBkgDstar));
507
508 for (int iEvt = 0; iEvt < 5; iEvt++) {
509 if (TMath::Abs(sPlotDatasetDstar->GetSWeight(iEvt, "nSignalDstar") + sPlotDatasetDstar->GetSWeight(iEvt, "nBkgDstar") - 1) > 5.e-3)
510 B2FATAL("Dstar: sPlot error: sum of weights not equal to 1");
511 }
512
513 RooDataSet* DstarDatasetSWeighted = new RooDataSet(DstarDataset->GetName(), DstarDataset->GetTitle(), DstarDataset,
514 *DstarDataset->get());
515
516 TTree* treeDstarSWeighted = DstarDatasetSWeighted->GetClonedTree();
517 treeDstarSWeighted->SetName("treeDstarSWeighted");
518
519 B2INFO("Dstar: sPlot done. Proceed to histogramming");
520 return treeDstarSWeighted;
521}
522
523std::unique_ptr<TList> SVDdEdxCalibrationAlgorithm::DstarHistogramming(TTree* inputTree)
524{
525 gROOT->SetBatch(true);
526 inputTree->SetEstimate(-1);
527 std::vector<double> pbins = CreatePBinningScheme();
528
529 TH2F* hDstarKMomentum = new TH2F("hist_d1_321_truncMomentum", "hist_d1_321_trunc;Momentum [GeV/c];dEdx [arb. units]", m_numPBins,
530 pbins.data(),
531 m_numDEdxBins, 0, m_dedxCutoff);
532 // the pion payload
533 TH2F* hDstarPiMomentum = new TH2F("hist_d1_211_truncMomentum", "hist_d1_211_trunc;Momentum [GeV/c];dEdx [arb. units]", m_numPBins,
534 pbins.data(),
535 m_numDEdxBins, 0, m_dedxCutoff);
536
537 inputTree->Draw("KaonSVDdEdx:KaonSVDdEdxTrackMomentum>>hist_d1_321_truncMomentum", "nSignalDstar_sw * (KaonSVDdEdx>0)", "goff");
538 // the pion one will be built from both pions in the Dstar decay tree
539 TH2F* hDstarPiPart1Momentum = (TH2F*)hDstarPiMomentum->Clone("hist_d1_211_truncPart1Momentum");
540 TH2F* hDstarPiPart2Momentum = (TH2F*)hDstarPiMomentum->Clone("hist_d1_211_truncPart2Momentum");
541
542 inputTree->Draw("PionDSVDdEdx:PionDSVDdEdxTrackMomentum>>hist_d1_211_truncPart1Momentum", "nSignalDstar_sw * (PionDSVDdEdx>0)",
543 "goff");
544 inputTree->Draw("SlowPionSVDdEdx:SlowPionSVDdEdxTrackMomentum>>hist_d1_211_truncPart2Momentum",
545 "nSignalDstar_sw * (SlowPionSVDdEdx>0)",
546 "goff");
547 hDstarPiMomentum->Add(hDstarPiPart1Momentum);
548 hDstarPiMomentum->Add(hDstarPiPart2Momentum);
549
550 inputTree->Draw(Form("KaonSVDdEdxTrackMomentum/%f", m_KaonPDGMass), "", "goff",
551 ((inputTree->GetEntries()) / m_numBGBins)*m_numBGBins);
552 double* KaonMomentumDataset = inputTree->GetV1();
553 TKDTreeBinning* kdBinsK = new TKDTreeBinning(((inputTree->GetEntries()) / m_numBGBins)*m_numBGBins, 1, KaonMomentumDataset,
554 m_numBGBins);
555 const double* binsMinEdgesKOriginal = kdBinsK->SortOneDimBinEdges();
556 double* binsMinEdgesK = const_cast<double*>(binsMinEdgesKOriginal);
557 binsMinEdgesK[0] = 0.1;
558 binsMinEdgesK[m_numBGBins + 1] = 50.;
559
560// get a distribution that contains both pions to get a typical kinematics for the binning scheme
561 inputTree->Draw(Form("SlowPionSVDdEdxTrackMomentum/%f * (event %% 2 == 0) + PionDSVDdEdxTrackMomentum/%f * (event %% 2 ==1)",
562 m_PionPDGMass, m_PionPDGMass), "", "goff",
563 ((inputTree->GetEntries()) / m_numBGBins)*m_numBGBins);
564 double* PionMomentumDataset = inputTree->GetV1();
565
566 TKDTreeBinning* kdBinsPi = new TKDTreeBinning(((inputTree->GetEntries()) / m_numBGBins)*m_numBGBins, 1, PionMomentumDataset,
567 m_numBGBins);
568 const double* binsMinEdgesPiOriginal = kdBinsPi->SortOneDimBinEdges();
569 double* binsMinEdgesPi = const_cast<double*>(binsMinEdgesPiOriginal);
570 binsMinEdgesPi[0] = 0.1;
571 binsMinEdgesPi[m_numBGBins + 1] = 50.;
572
573 TH2F* hDstarKBetaGamma = new TH2F("hist_d1_321_truncBetaGamma", "hist_d1_321_truncBetaGamma;#beta*#gamma;dEdx [arb. units]",
574 m_numBGBins,
575 binsMinEdgesK,
576 m_numDEdxBins, 0, m_dedxMaxPossible);
577 // the pion payload
578 TH2F* hDstarPiBetaGamma = new TH2F("hist_d1_211_truncBetaGamma", "hist_d1_211_truncBetaGamma;#beta*#gamma;dEdx [arb. units]",
579 m_numBGBins,
580 binsMinEdgesPi,
581 m_numDEdxBins, 0, m_dedxMaxPossible);
582
583 inputTree->Draw(Form("KaonSVDdEdx:KaonSVDdEdxTrackMomentum/%f>>hist_d1_321_truncBetaGamma", m_KaonPDGMass),
584 "nSignalDstar_sw * (KaonSVDdEdx>0)", "goff");
585 // the pion one will be built from both pions in the Dstar decay tree
586 TH2F* hDstarPiPart1BetaGamma = (TH2F*)hDstarPiBetaGamma->Clone("hist_d1_211_truncPart1BetaGamma");
587 TH2F* hDstarPiPart2BetaGamma = (TH2F*)hDstarPiBetaGamma->Clone("hist_d1_211_truncPart2BetaGamma");
588
589 inputTree->Draw(Form("PionDSVDdEdx:PionDSVDdEdxTrackMomentum/%f>>hist_d1_211_truncPart1BetaGamma", m_PionPDGMass),
590 "nSignalDstar_sw * (PionDSVDdEdx>0)",
591 "goff");
592 inputTree->Draw(Form("SlowPionSVDdEdx:SlowPionSVDdEdxTrackMomentum/%f>>hist_d1_211_truncPart2BetaGamma", m_PionPDGMass),
593 "nSignalDstar_sw * (SlowPionSVDdEdx>0)", "goff");
594 hDstarPiBetaGamma->Add(hDstarPiPart1BetaGamma);
595 hDstarPiBetaGamma->Add(hDstarPiPart2BetaGamma);
596
597
598
599 // produce the 1D profiles
600
601
602 TH1F* PionProfileMomentum = (TH1F*)hDstarPiMomentum->ProfileX("PionProfileMomentum");
603 PionProfileMomentum->SetTitle("PionProfile");
604 PionProfileMomentum->GetYaxis()->SetRangeUser(0, m_dedxCutoff);
605 PionProfileMomentum->GetXaxis()->SetTitle("Momentum, GeV/c");
606 PionProfileMomentum->GetYaxis()->SetTitle("dE/dx");
607 PionProfileMomentum->SetLineColor(kRed);
608
609 TH1F* PionProfileBetaGamma = (TH1F*)hDstarPiBetaGamma->ProfileX("PionProfileBetaGamma");
610 if (m_CustomProfile) {
611 PionProfileBetaGamma = PrepareProfile(hDstarPiBetaGamma, "PionProfileBetaGamma");
612 }
613 PionProfileBetaGamma->SetTitle("PionProfile");
614 PionProfileBetaGamma->GetYaxis()->SetRangeUser(0, m_dedxCutoff);
615 PionProfileBetaGamma->GetXaxis()->SetTitle("#beta*#gamma");
616 PionProfileBetaGamma->GetYaxis()->SetTitle("dE/dx");
617 PionProfileBetaGamma->SetLineColor(kRed);
618
619
620 TH1F* KaonProfileMomentum = (TH1F*)hDstarKMomentum->ProfileX("KaonProfileMomentum");
621 KaonProfileMomentum->SetTitle("KaonProfile");
622 KaonProfileMomentum->GetYaxis()->SetRangeUser(0, m_dedxCutoff);
623 KaonProfileMomentum->GetXaxis()->SetTitle("Momentum, GeV/c");
624 KaonProfileMomentum->GetYaxis()->SetTitle("dE/dx");
625 KaonProfileMomentum->SetLineColor(kRed);
626
627
628 TH1F* KaonProfileBetaGamma = (TH1F*)hDstarKBetaGamma->ProfileX("KaonProfileBetaGamma");
629 if (m_CustomProfile) {
630 KaonProfileBetaGamma = PrepareProfile(hDstarKBetaGamma, "KaonProfileBetaGamma");
631 }
632 KaonProfileBetaGamma->SetTitle("KaonProfile");
633
634 KaonProfileBetaGamma->GetYaxis()->SetRangeUser(0, m_dedxCutoff);
635 KaonProfileBetaGamma->GetXaxis()->SetTitle("#beta*#gamma");
636 KaonProfileBetaGamma->GetYaxis()->SetTitle("dE/dx");
637 KaonProfileBetaGamma->SetLineColor(kRed);
638
639
640
641 // normalisation
642 hDstarKMomentum->Sumw2();
643 hDstarKMomentum = Normalise2DHisto(hDstarKMomentum);
644
645 hDstarPiMomentum->Sumw2();
646 hDstarPiMomentum = Normalise2DHisto(hDstarPiMomentum);
647
648 std::unique_ptr<TList> histList(new TList);
649 histList->Add(KaonProfileMomentum);
650 histList->Add(KaonProfileBetaGamma);
651 histList->Add(hDstarKMomentum);
652
653 histList->Add(PionProfileMomentum);
654 histList->Add(PionProfileBetaGamma);
655 histList->Add(hDstarPiMomentum);
656
657 if (m_isMakePlots) {
658 TFile DstarHistogrammingPlotFile("SVDdEdxCalibrationDstarHistogramming.root", "RECREATE");
659 histList->Write();
660 DstarHistogrammingPlotFile.Close();
661 }
662
663 return histList;
664}
665
666
667std::unique_ptr<TList> SVDdEdxCalibrationAlgorithm::GammaHistogramming(std::shared_ptr<TTree> preselTree)
668{
669 B2INFO("Histogramming the converted photon selection...");
670 gROOT->SetBatch(true);
671
672
673 if (preselTree->GetEntries() == 0) {
674 B2FATAL("The Gamma tree is empty, stopping here");
675 }
676 preselTree->SetEstimate(-1);
677 std::vector<double> pbins = CreatePBinningScheme();
678
679
680 TH2F* hGammaEMomentum = new TH2F("hist_d1_11_truncMomentum", "hist_d1_11_trunc;Momentum [GeV/c];dEdx [arb. units]", m_numPBins,
681 pbins.data(), m_numDEdxBins, 0, m_dedxCutoff);
682
683 TH2F* hGammaEPart1Momentum = (TH2F*)hGammaEMomentum->Clone("hist_d1_11_truncPart1Momentum");
684 TH2F* hGammaEPart2Momentum = (TH2F*)hGammaEMomentum->Clone("hist_d1_11_truncPart2Momentum");
685
686 preselTree->Draw("FirstElectronSVDdEdx:FirstElectronSVDdEdxTrackMomentum>>hist_d1_11_truncPart1Momentum",
687 "FirstElectronSVDdEdx>0 && DIRA>0.995 && dr>1.2", "goff");
688 preselTree->Draw("SecondElectronSVDdEdx:SecondElectronSVDdEdxTrackMomentum>>hist_d1_11_truncPart2Momentum",
689 "SecondElectronSVDdEdx>0 && DIRA>0.995 && dr>1.2", "goff");
690 hGammaEMomentum->Add(hGammaEPart1Momentum);
691 hGammaEMomentum->Add(hGammaEPart2Momentum);
692
693// get a distribution that contains both pions to get a typical kinematics for the binning scheme
694 preselTree->Draw(
695 Form("FirstElectronSVDdEdxTrackMomentum/%f* (event %% 2==0) + SecondElectronSVDdEdxTrackMomentum/%f* (event %% 2==1)",
696 m_ElectronPDGMass, m_ElectronPDGMass),
697 "", "goff", ((preselTree->GetEntries()) / m_numBGBins)*m_numBGBins);
698 double* ElectronMomentumDataset = preselTree->GetV1();
699
700 TKDTreeBinning* kdBinsE = new TKDTreeBinning(((preselTree->GetEntries()) / m_numBGBins)*m_numBGBins, 1, ElectronMomentumDataset,
701 m_numBGBins);
702 const double* binsMinEdgesEOriginal = kdBinsE->SortOneDimBinEdges();
703 double* binsMinEdgesE = const_cast<double*>(binsMinEdgesEOriginal);
704 binsMinEdgesE[0] = 0.;
705 binsMinEdgesE[m_numBGBins + 1] = 10000.;
706
707
708 TH2F* hGammaEBetaGamma = new TH2F("hist_d1_11_truncBetaGamma", "hist_d1_11_truncBetaGamma;#beta*#gamma;dEdx [arb. units]",
709 m_numBGBins,
710 binsMinEdgesE,
711 m_numDEdxBins, 0, m_dedxMaxPossible);
712 TH2F* hGammaEPart1BetaGamma = (TH2F*)hGammaEBetaGamma->Clone("hist_d1_11_truncPart1BetaGamma");
713 TH2F* hGammaEPart2BetaGamma = (TH2F*)hGammaEBetaGamma->Clone("hist_d1_11_truncPart2BetaGamma");
714
715 preselTree->Draw(Form("FirstElectronSVDdEdx:FirstElectronSVDdEdxTrackMomentum/%f>>hist_d1_11_truncPart1BetaGamma",
716 m_ElectronPDGMass),
717 "FirstElectronSVDdEdx>0 && DIRA>0.995 && dr>1.2 && FirstElectronSVDdEdx<1.8e6", "goff");
718 preselTree->Draw(Form("SecondElectronSVDdEdx:SecondElectronSVDdEdxTrackMomentum/%f>>hist_d1_11_truncPart2BetaGamma",
719 m_ElectronPDGMass),
720 "SecondElectronSVDdEdx>0 && DIRA>0.995 && dr>1.2 && SecondElectronSVDdEdx<1.8e6", "goff");
721 hGammaEBetaGamma->Add(hGammaEPart1BetaGamma);
722 hGammaEBetaGamma->Add(hGammaEPart2BetaGamma);
723
724
725 // produce the 1D profile (for data-MC comparisons)
726 TH1F* ElectronProfileMomentum = (TH1F*)hGammaEMomentum->ProfileX("ElectronProfileMomentum");
727 ElectronProfileMomentum->SetTitle("ElectronProfile");
728 ElectronProfileMomentum->GetYaxis()->SetRangeUser(0, m_dedxCutoff);
729 ElectronProfileMomentum->GetXaxis()->SetTitle("Momentum, GeV/c");
730 ElectronProfileMomentum->GetYaxis()->SetTitle("dE/dx");
731 ElectronProfileMomentum->SetLineColor(kRed);
732
733// beta*gamma profile for the fit
734 TH1F* ElectronProfileBetaGamma = (TH1F*)hGammaEBetaGamma->ProfileX("ElectronProfileBetaGamma");
735 if (m_CustomProfile) {
736 ElectronProfileBetaGamma = PrepareProfile(hGammaEBetaGamma, "ElectronProfileBetaGamma");
737 }
738 ElectronProfileBetaGamma->SetTitle("ElectronProfile");
739
740 ElectronProfileBetaGamma->GetYaxis()->SetRangeUser(0, m_dedxCutoff);
741 ElectronProfileBetaGamma->GetXaxis()->SetTitle("#beta*#gamma");
742 ElectronProfileBetaGamma->GetYaxis()->SetTitle("dE/dx");
743 ElectronProfileBetaGamma->SetLineColor(kRed);
744
745
746 hGammaEMomentum = Normalise2DHisto(hGammaEMomentum);
747
748 std::unique_ptr<TList> histList(new TList);
749 histList->Add(ElectronProfileMomentum);
750 histList->Add(ElectronProfileBetaGamma);
751 histList->Add(hGammaEMomentum);
752
753 if (m_isMakePlots) {
754 TFile GammaHistogrammingPlotFile("SVDdEdxCalibrationGammaHistogramming.root", "RECREATE");
755 histList->Write();
756 GammaHistogrammingPlotFile.Close();
757 }
758
759 return histList;
760
761}
762
763std::unique_ptr<TList> SVDdEdxCalibrationAlgorithm::GenerateNewHistograms(std::shared_ptr<TTree> ttreeLambda,
764 std::shared_ptr<TTree> ttreeDstar,
765 std::shared_ptr<TTree> ttreeGamma, std::shared_ptr<TTree> ttreeGeneric)
766{
767 gROOT->SetBatch(true);
768 gStyle->SetOptStat(0);
769// run the background subtraction and histogramming parts
770 TTree* treeLambda = LambdaMassFit(ttreeLambda);
771 std::unique_ptr<TList> HistListLambda = LambdaHistogramming(treeLambda);
772 TH1F* ProtonProfileBetaGamma = (TH1F*) HistListLambda->FindObject("ProtonProfileBetaGamma");
773 TH2F* Proton2DHistogram = (TH2F*) HistListLambda->FindObject("hist_d1_2212_truncMomentum");
774
775 TTree* treeDstar = DstarMassFit(ttreeDstar);
776 std::unique_ptr<TList> HistListDstar = DstarHistogramming(treeDstar);
777 TH1F* PionProfileBetaGamma = (TH1F*) HistListDstar->FindObject("PionProfileBetaGamma");
778 TH2F* Pion2DHistogram = (TH2F*) HistListDstar->FindObject("hist_d1_211_truncMomentum");
779 TH1F* KaonProfileBetaGamma = (TH1F*) HistListDstar->FindObject("KaonProfileBetaGamma");
780 TH2F* Kaon2DHistogram = (TH2F*) HistListDstar->FindObject("hist_d1_321_truncMomentum");
781
782 std::unique_ptr<TList> HistListGamma = GammaHistogramming(ttreeGamma);
783 TH1F* ElectronProfileBetaGamma = (TH1F*) HistListGamma->FindObject("ElectronProfileBetaGamma");
784 TH2F* Electron2DHistogram = (TH2F*) HistListGamma->FindObject("hist_d1_11_truncMomentum");
785
786 int cred = TColor::GetColor("#e31a1c");
787 PionProfileBetaGamma->SetMarkerSize(4);
788 PionProfileBetaGamma->SetLineWidth(2);
789 PionProfileBetaGamma->SetMarkerColor(cred);
790 PionProfileBetaGamma->SetLineColor(cred);
791
792 int cpink = TColor::GetColor("#807dba");
793 KaonProfileBetaGamma->SetMarkerSize(4);
794 KaonProfileBetaGamma->SetLineWidth(2);
795 KaonProfileBetaGamma->SetMarkerColor(cpink);
796 KaonProfileBetaGamma->SetLineColor(cpink);
797
798 int cblue = TColor::GetColor("#084594");
799 ProtonProfileBetaGamma->SetMarkerSize(4);
800 ProtonProfileBetaGamma->SetLineWidth(2);
801 ProtonProfileBetaGamma->SetMarkerColor(cblue);
802 ProtonProfileBetaGamma->SetLineColor(cblue);
803
804 int cgreen = TColor::GetColor("#238b45");
805 ElectronProfileBetaGamma->SetMarkerSize(4);
806 ElectronProfileBetaGamma->SetLineWidth(2);
807 ElectronProfileBetaGamma->SetMarkerColor(cgreen);
808 ElectronProfileBetaGamma->SetLineColor(cgreen);
809
810// prepare the fitting
811
812 PionProfileBetaGamma->GetYaxis()->SetRangeUser(5.e5, 5.5e6);
813 KaonProfileBetaGamma->GetYaxis()->SetRangeUser(5.e5, 5.5e6);
814 ProtonProfileBetaGamma->GetYaxis()->SetRangeUser(5.e5, 5.5e6);
815
816// enhance the proton histogram (which ends at beta*gamma around 3) by adding pion data above this value – otherwise the fit is very unstable
817 auto PionEdges = PionProfileBetaGamma->GetXaxis()->GetXbins()->GetArray();
818 auto ProtonEdges = ProtonProfileBetaGamma->GetXaxis()->GetXbins()->GetArray();
819
820 std::vector<float> CombinedEdgesVector;
821
822 double borderline = 3.;
823
824 for (int i = 0; i < ProtonProfileBetaGamma->GetNbinsX() + 1; i++)
825 if (ProtonEdges[i] < borderline) CombinedEdgesVector.push_back(ProtonEdges[i]);
826
827
828 for (int i = 0; i < PionProfileBetaGamma->GetNbinsX() + 1; i++)
829 if (PionEdges[i] > borderline) CombinedEdgesVector.push_back(PionEdges[i]);
830
831
832 TH1F* CombinedHistogramPAndPi = new TH1F("CombinedHistogramPAndPi", "histo_for_fit", CombinedEdgesVector.size() - 1,
833 CombinedEdgesVector.data());
834
835 int iterator = 1;
836 for (int i = 1; i < ProtonProfileBetaGamma->GetNbinsX() + 1; i++)
837 if (ProtonEdges[i - 1] < borderline) {
838 CombinedHistogramPAndPi->SetBinContent(i, ProtonProfileBetaGamma->GetBinContent(i));
839 CombinedHistogramPAndPi->SetBinError(i, ProtonProfileBetaGamma->GetBinError(i));
840 iterator++;
841 }
842
843 for (int i = 1; i < PionProfileBetaGamma->GetNbinsX() + 1; i++)
844 if (PionEdges[i - 1] > borderline) {
845
846 CombinedHistogramPAndPi->SetBinContent(iterator, PionProfileBetaGamma->GetBinContent(i));
847 CombinedHistogramPAndPi->SetBinError(iterator, PionProfileBetaGamma->GetBinError(i));
848 iterator++;
849 }
850
851// define the beta*gamma vs momentum function
852 TF1* BetaGammaFunctionPion = new TF1("BetaGammaFunctionPion", "[0] + [1] * x/[2] + [5]/(x^2/[2]^2 + [3])**[4] + [6]* (x/[2])**0.5",
853 0.01, 25.);
854
855 BetaGammaFunctionPion->SetNpx(1000);
856
857 BetaGammaFunctionPion->SetParameters(5.e5, 2.e3, 1, 0.15, 1.2, 6.e5, 3.e5);
858
859 BetaGammaFunctionPion->SetParLimits(0, 3.e5, 7.e5);
860 BetaGammaFunctionPion->SetParLimits(1, -3.e4, 1.e4);
861 BetaGammaFunctionPion->SetParLimits(3, 0.1, 0.2);
862 BetaGammaFunctionPion->SetParLimits(4, 0.9, 1.6);
863 BetaGammaFunctionPion->SetParLimits(5, 3.e5, 7.e5);
864 BetaGammaFunctionPion->SetParLimits(6, 0., 1.e6);
865 BetaGammaFunctionPion->FixParameter(2, 1);
866 if (m_FixUnstableFitParameter) BetaGammaFunctionPion->FixParameter(3, 0.15);
867
868// fit it to the pion data
869 ROOT::Math::MinimizerOptions::SetDefaultMinimizer("Minuit2", "Migrad");
870 auto FitResultBetaGammaPion = PionProfileBetaGamma->Fit("BetaGammaFunctionPion", "0SI", "", 0.4, 25);
871
872 if ((FitResultBetaGammaPion->Status() > 1) || (BetaGammaFunctionPion->Eval(1) < 5.e5) || (BetaGammaFunctionPion->Eval(1) > 5.e6)) {
873 BetaGammaFunctionPion->FixParameter(3, 0.15);
874 FitResultBetaGammaPion = PionProfileBetaGamma->Fit("BetaGammaFunctionPion", "0SI", "", 0.4, 25);
875 }
876 if ((FitResultBetaGammaPion->Status() > 1) || (BetaGammaFunctionPion->Eval(1) < 5.e5) || (BetaGammaFunctionPion->Eval(1) > 5.e6)) {
877 BetaGammaFunctionPion->FixParameter(3, 0.15);
878 FitResultBetaGammaPion = PionProfileBetaGamma->Fit("BetaGammaFunctionPion", "0SI", "", 0.45, 25);
879 }
880 if ((FitResultBetaGammaPion->Status() > 1) || (BetaGammaFunctionPion->Eval(1) < 5.e5) || (BetaGammaFunctionPion->Eval(1) > 5.e6)) {
881 BetaGammaFunctionPion->FixParameter(3, 0.15);
882 FitResultBetaGammaPion = PionProfileBetaGamma->Fit("BetaGammaFunctionPion", "0S", "", 0.5, 25);
883 }
884
885 B2INFO("BetaGamma fit for pions done. Fit status: " << FitResultBetaGammaPion->Status());
886 // B2INFO(FitResultBetaGammaPion->Print(Belle2::LogConfig::c_Info));
887 B2INFO("Fit parameters:");
888 B2INFO("p0: " << BetaGammaFunctionPion->GetParameter(0) << " +- " << BetaGammaFunctionPion->GetParError(0));
889 B2INFO("p1: " << BetaGammaFunctionPion->GetParameter(1) << " +- " << BetaGammaFunctionPion->GetParError(1));
890 B2INFO("p2: " << BetaGammaFunctionPion->GetParameter(2) << " +- " << BetaGammaFunctionPion->GetParError(2));
891 B2INFO("p3: " << BetaGammaFunctionPion->GetParameter(3) << " +- " << BetaGammaFunctionPion->GetParError(3));
892 B2INFO("p4: " << BetaGammaFunctionPion->GetParameter(4) << " +- " << BetaGammaFunctionPion->GetParError(4));
893 B2INFO("p5: " << BetaGammaFunctionPion->GetParameter(5) << " +- " << BetaGammaFunctionPion->GetParError(5));
894 B2INFO("p6: " << BetaGammaFunctionPion->GetParameter(6) << " +- " << BetaGammaFunctionPion->GetParError(6));
895
896 // repeat the same for kaons
897 TF1* BetaGammaFunctionKaon = new TF1("BetaGammaFunctionKaon", "[0] + [1] * x/[2] + [5]/(x^2/[2]^2 + [3])**[4]+ [6]* (x/[2])**0.5",
898 0.01, 25.);
899
900 BetaGammaFunctionKaon->SetNpx(1000);
901 BetaGammaFunctionKaon->SetParameters(5.e5, 2.e3, 1, 0.15, 1.2, 6.e5, 3.e5);
902
903 BetaGammaFunctionKaon->SetParLimits(0, 3.e5, 7.e5);
904 BetaGammaFunctionKaon->SetParLimits(1, -3.e4, 1.e4);
905 BetaGammaFunctionKaon->SetParLimits(3, 0.1, 0.2);
906 BetaGammaFunctionKaon->SetParLimits(4, 0.9, 1.6);
907 BetaGammaFunctionKaon->SetParLimits(5, 3.e5, 7.e5);
908 BetaGammaFunctionKaon->SetParLimits(6, 0., 1.e6);
909
910 BetaGammaFunctionKaon->FixParameter(2, 1);
911 if (m_FixUnstableFitParameter) BetaGammaFunctionKaon->FixParameter(3, 0.15);
912
913 BetaGammaFunctionKaon->SetLineColor(KaonProfileBetaGamma->GetMarkerColor());
914
915 auto FitResultBetaGammaKaon = KaonProfileBetaGamma->Fit("BetaGammaFunctionKaon", "0SI", "", 0.4, 8.5);
916
917 if ((FitResultBetaGammaKaon->Status() > 1) || (BetaGammaFunctionKaon->Eval(1) < 5.e5) || (BetaGammaFunctionKaon->Eval(1) > 5.e6)) {
918 BetaGammaFunctionKaon->FixParameter(3, 0.15);
919 FitResultBetaGammaKaon = KaonProfileBetaGamma->Fit("BetaGammaFunctionKaon", "0SI", "", 0.4, 8.5);
920 }
921 if ((FitResultBetaGammaKaon->Status() > 1) || (BetaGammaFunctionKaon->Eval(1) < 5.e5) || (BetaGammaFunctionKaon->Eval(1) > 5.e6)) {
922 BetaGammaFunctionKaon->FixParameter(3, 0.15);
923 FitResultBetaGammaKaon = KaonProfileBetaGamma->Fit("BetaGammaFunctionKaon", "0SI", "", 0.45, 8);
924 }
925 if ((FitResultBetaGammaKaon->Status() > 1) || (BetaGammaFunctionKaon->Eval(1) < 5.e5) || (BetaGammaFunctionKaon->Eval(1) > 5.e6)) {
926 BetaGammaFunctionKaon->FixParameter(3, 0.15);
927 FitResultBetaGammaKaon = KaonProfileBetaGamma->Fit("BetaGammaFunctionKaon", "0S", "", 0.5, 8);
928 }
929
930 B2INFO("BetaGamma fit for kaons done. Fit status: " << FitResultBetaGammaKaon->Status());
931 B2INFO("Fit parameters:");
932 B2INFO("p0: " << BetaGammaFunctionKaon->GetParameter(0) << " +- " << BetaGammaFunctionKaon->GetParError(0));
933 B2INFO("p1: " << BetaGammaFunctionKaon->GetParameter(1) << " +- " << BetaGammaFunctionKaon->GetParError(1));
934 B2INFO("p2: " << BetaGammaFunctionKaon->GetParameter(2) << " +- " << BetaGammaFunctionKaon->GetParError(2));
935 B2INFO("p3: " << BetaGammaFunctionKaon->GetParameter(3) << " +- " << BetaGammaFunctionKaon->GetParError(3));
936 B2INFO("p4: " << BetaGammaFunctionKaon->GetParameter(4) << " +- " << BetaGammaFunctionKaon->GetParError(4));
937 B2INFO("p5: " << BetaGammaFunctionKaon->GetParameter(5) << " +- " << BetaGammaFunctionKaon->GetParError(5));
938 B2INFO("p6: " << BetaGammaFunctionKaon->GetParameter(6) << " +- " << BetaGammaFunctionKaon->GetParError(6));
939
940 // repeat the same for protons
941 TF1* BetaGammaFunctionProton = new TF1("BetaGammaFunctionProton",
942 "[0] + [1] * x/[2] + [5]/(x^2/[2]^2 + [3])**[4]+ [6]* (x/[2])**0.5", 0.01, 25.);
943
944 BetaGammaFunctionProton->SetNpx(1000);
945
946 BetaGammaFunctionProton->SetParameters(5.e5, 2.e3, 1, 0.15, 1.2, 6.e5, 3.e5);
947
948 BetaGammaFunctionProton->SetParLimits(0, 3.e5, 7.e5);
949 BetaGammaFunctionProton->SetParLimits(1, -3.e4, 1.e4);
950 BetaGammaFunctionProton->SetParLimits(3, 0.1, 0.2);
951 BetaGammaFunctionProton->SetParLimits(4, 0.9, 1.6);
952 BetaGammaFunctionProton->SetParLimits(5, 3.e5, 7.e5);
953 BetaGammaFunctionProton->SetParLimits(6, 0., 1.e6);
954
955 BetaGammaFunctionProton->FixParameter(2, 1);
956 if (m_FixUnstableFitParameter) BetaGammaFunctionProton->FixParameter(3, 0.15);
957
958 BetaGammaFunctionProton->SetLineColor(ProtonProfileBetaGamma->GetMarkerColor());
959
960 auto FitResultBetaGammaProton = CombinedHistogramPAndPi->Fit("BetaGammaFunctionProton", "0SI", "", 0.45, 15);
961
962 if ((FitResultBetaGammaProton->Status() > 1) || (BetaGammaFunctionProton->Eval(1) < 5.e5)
963 || (BetaGammaFunctionProton->Eval(1) > 5.e6)) {
964 BetaGammaFunctionProton->FixParameter(3, 0.15);
965 FitResultBetaGammaProton = CombinedHistogramPAndPi->Fit("BetaGammaFunctionProton", "0SI", "", 0.45, 15);
966 }
967 if ((FitResultBetaGammaProton->Status() > 1) || (BetaGammaFunctionProton->Eval(1) < 5.e5)
968 || (BetaGammaFunctionProton->Eval(1) > 5.e6)) {
969 BetaGammaFunctionProton->FixParameter(3, 0.15);
970 FitResultBetaGammaProton = CombinedHistogramPAndPi->Fit("BetaGammaFunctionProton", "0SI", "", 0.45, 10);
971 }
972 if ((FitResultBetaGammaProton->Status() > 1) || (BetaGammaFunctionProton->Eval(1) < 5.e5)
973 || (BetaGammaFunctionProton->Eval(1) > 5.e6)) {
974 BetaGammaFunctionProton->FixParameter(3, 0.15);
975 FitResultBetaGammaProton = CombinedHistogramPAndPi->Fit("BetaGammaFunctionProton", "0S", "", 0.5, 10);
976 }
977
978 B2INFO("BetaGamma fit for protons done. Fit status: " << FitResultBetaGammaProton->Status());
979 B2INFO("Fit parameters:");
980 B2INFO("p0: " << BetaGammaFunctionProton->GetParameter(0) << " +- " << BetaGammaFunctionProton->GetParError(0));
981 B2INFO("p1: " << BetaGammaFunctionProton->GetParameter(1) << " +- " << BetaGammaFunctionProton->GetParError(1));
982 B2INFO("p2: " << BetaGammaFunctionProton->GetParameter(2) << " +- " << BetaGammaFunctionProton->GetParError(2));
983 B2INFO("p3: " << BetaGammaFunctionProton->GetParameter(3) << " +- " << BetaGammaFunctionProton->GetParError(3));
984 B2INFO("p4: " << BetaGammaFunctionProton->GetParameter(4) << " +- " << BetaGammaFunctionProton->GetParError(4));
985 B2INFO("p5: " << BetaGammaFunctionProton->GetParameter(5) << " +- " << BetaGammaFunctionProton->GetParError(5));
986 B2INFO("p6: " << BetaGammaFunctionProton->GetParameter(6) << " +- " << BetaGammaFunctionProton->GetParError(6));
987
988 if (m_isMakePlots) {
989// plot a summary of all beta*gamma fits for hadrons
990 std::unique_ptr<TCanvas> CombinedCanvasHadrons(new TCanvas("CombinedCanvasHadrons", "Hadron beta*gamma fits", 10, 10, 1000, 700));
991 gStyle->SetOptFit(1111);
992
993 PionProfileBetaGamma->Draw();
994 PionProfileBetaGamma->GetListOfFunctions()->Add(BetaGammaFunctionPion);
995 KaonProfileBetaGamma->Draw("SAME");
996 KaonProfileBetaGamma->GetListOfFunctions()->Add(BetaGammaFunctionKaon);
997 ProtonProfileBetaGamma->Draw("SAME");
998 ProtonProfileBetaGamma->GetListOfFunctions()->Add(BetaGammaFunctionProton);
999// BetaGammaFunctionPion->Draw("SAME");
1000// BetaGammaFunctionKaon->Draw("SAME");
1001// BetaGammaFunctionProton->Draw("SAME");
1002 auto legend = new TLegend(0.4, 0.7, 0.65, 0.9);
1003 legend->AddEntry(PionProfileBetaGamma, "Pions", "lep");
1004 legend->AddEntry(KaonProfileBetaGamma, "Kaons", "lep");
1005 legend->AddEntry(ProtonProfileBetaGamma, "Protons", "lep");
1006 legend->Draw();
1007
1008 gPad->SetLogx();
1009 gPad->SetLogy();
1010
1011 CombinedCanvasHadrons->Print("HadronBetaGammaFits.pdf");
1012 TFile HadronFitPlotFile("SVDdEdxCalibrationHadronFitPlotFile.root", "RECREATE");
1013 PionProfileBetaGamma->Write();
1014 KaonProfileBetaGamma->Write();
1015 ProtonProfileBetaGamma->Write();
1016 CombinedHistogramPAndPi->Write();
1017 BetaGammaFunctionPion->Write();
1018 BetaGammaFunctionKaon->Write();
1019 BetaGammaFunctionProton->Write();
1020 CombinedCanvasHadrons->Write();
1021 HadronFitPlotFile.Close();
1022 }
1023
1024
1025 // in case we assume that all hadrons are equal
1026 if (m_UsePionBGFunctionForEverything) {
1027 BetaGammaFunctionKaon = (TF1*) BetaGammaFunctionPion->Clone("BetaGammaFunctionKaon");
1028 BetaGammaFunctionProton = (TF1*) BetaGammaFunctionPion->Clone("BetaGammaFunctionProton");
1029 }
1030 if (m_UseProtonBGFunctionForEverything) {
1031 BetaGammaFunctionKaon = (TF1*) BetaGammaFunctionProton->Clone("BetaGammaFunctionKaon");
1032 BetaGammaFunctionPion = (TF1*) BetaGammaFunctionProton->Clone("BetaGammaFunctionPion");
1033 }
1034
1035 // sanity checks: are all fits ok?
1036 if ((FitResultBetaGammaProton->Status() > 1) || (BetaGammaFunctionProton->Eval(1) < 5.e5)
1037 || (BetaGammaFunctionProton->Eval(1) > 5.e6)) {
1038 if (FitResultBetaGammaPion->Status() == 0) {
1039 BetaGammaFunctionProton = (TF1*) BetaGammaFunctionPion->Clone("BetaGammaFunctionProton");
1040 } else if (FitResultBetaGammaKaon->Status() == 0) {
1041 BetaGammaFunctionProton = (TF1*) BetaGammaFunctionKaon->Clone("BetaGammaFunctionProton");
1042 } else {
1043 B2WARNING("Problem with the beta*gamma fit for protons, reverting to the default values");
1044 BetaGammaFunctionProton->SetParameters(450258, -10900.8, 1, 0.126797, 1.155, 641907, 86304.5);
1045 }
1046 }
1047
1048 if ((FitResultBetaGammaKaon->Status() > 1) || (BetaGammaFunctionKaon->Eval(1) < 5.e5) || (BetaGammaFunctionKaon->Eval(1) > 5.e6)) {
1049 if (FitResultBetaGammaProton->Status() == 0) {
1050 BetaGammaFunctionKaon = (TF1*) BetaGammaFunctionProton->Clone("BetaGammaFunctionKaon");
1051 } else if (FitResultBetaGammaPion->Status() == 0) {
1052 BetaGammaFunctionKaon = (TF1*) BetaGammaFunctionPion->Clone("BetaGammaFunctionKaon");
1053 } else {
1054 B2WARNING("Problem with the beta*gamma fit for kaons, reverting to the default values");
1055 BetaGammaFunctionKaon->SetParameters(543386, 3013.81, 1, 0.135517, 1.19742, 619509, 15484.4);
1056 }
1057 }
1058
1059 if ((FitResultBetaGammaPion->Status() > 1) || (BetaGammaFunctionPion->Eval(1) < 5.e5) || (BetaGammaFunctionPion->Eval(1) > 5.e6)) {
1060 if (FitResultBetaGammaKaon->Status() == 0) {
1061 BetaGammaFunctionPion = (TF1*) BetaGammaFunctionKaon->Clone("BetaGammaFunctionPion");
1062 } else if (FitResultBetaGammaProton->Status() == 0) {
1063 BetaGammaFunctionPion = (TF1*) BetaGammaFunctionProton->Clone("BetaGammaFunctionPion");
1064 } else {
1065 B2WARNING("Problem with the beta*gamma fit for pions, reverting to the default values");
1066 BetaGammaFunctionPion->SetParameters(537623, -1937.62, 1, 0.15292, 1.23803, 623678, 30400.9);
1067 }
1068 }
1069
1070
1071// electrons
1072
1073 TF1* BetaGammaFunctionElectron = new TF1("BetaGammaFunctionElectron", "[0] + [1]* x", 1, 10000.);
1074 BetaGammaFunctionElectron->SetParameters(6.e5, -1);
1075 BetaGammaFunctionElectron->SetParLimits(0, 3.e5, 8.e5);
1076 BetaGammaFunctionElectron->SetParLimits(1, -1.e5, 1.e5);
1077 auto FitResultBetaGammaElectron = ElectronProfileBetaGamma->Fit("BetaGammaFunctionElectron", "0SI", "", 100, 8000);
1078
1079
1080 if ((FitResultBetaGammaElectron->Status() > 1) || (BetaGammaFunctionElectron->Eval(1) < 3.e5)
1081 || (BetaGammaFunctionElectron->Eval(1) > 5.e6)) {
1082 FitResultBetaGammaElectron = ElectronProfileBetaGamma->Fit("BetaGammaFunctionElectron", "0S", "", 100, 10000);
1083 }
1084 B2INFO("BetaGamma fit for electrons done. Fit status: " << FitResultBetaGammaElectron->Status());
1085 B2INFO("Fit parameters:");
1086 B2INFO("p0: " << BetaGammaFunctionElectron->GetParameter(0) << " +- " << BetaGammaFunctionElectron->GetParError(0));
1087 B2INFO("p1: " << BetaGammaFunctionElectron->GetParameter(1) << " +- " << BetaGammaFunctionElectron->GetParError(1));
1088
1089
1090 ElectronProfileBetaGamma->SetMarkerSize(4);
1091 ElectronProfileBetaGamma->SetLineWidth(2);
1092 ElectronProfileBetaGamma->GetYaxis()->SetRangeUser(5e5, 1e6);
1093 ElectronProfileBetaGamma->GetListOfFunctions()->Add(BetaGammaFunctionElectron);
1094
1095 if (m_isMakePlots) {
1096 std::unique_ptr<TCanvas> ElectronCanvas(new TCanvas("ElectronCanvas", "Electron histogram", 10, 10, 1000, 700));
1097 ElectronProfileBetaGamma->Draw();
1098
1099 gPad->SetLogx();
1100
1101 ElectronCanvas->Print("ElectronBetaGammaFits.pdf");
1102 TFile ElectronFitPlotFile("SVDdEdxCalibrationElectronFitPlotFile.root", "RECREATE");
1103 ElectronProfileBetaGamma->Write();
1104 BetaGammaFunctionElectron->Write();
1105 ElectronCanvas->Write();
1106 ElectronFitPlotFile.Close();
1107 }
1108
1109 TF1* MomentumFunctionElectron = (TF1*) BetaGammaFunctionElectron->Clone("MomentumFunctionElectron");
1110 MomentumFunctionElectron->SetParameter(2, m_ElectronPDGMass);
1111 MomentumFunctionElectron->SetRange(0.01, 5.5);
1112 MomentumFunctionElectron->SetLineColor(kRed);
1113 MomentumFunctionElectron->SetLineWidth(4);
1114
1115 TF1* MomentumFunctionPion = (TF1*) BetaGammaFunctionPion->Clone("MomentumFunctionPion");
1116 MomentumFunctionPion->SetParameter(2, m_PionPDGMass);
1117 MomentumFunctionPion->SetRange(0.01, 5.5);
1118 MomentumFunctionPion->SetLineColor(kRed);
1119 MomentumFunctionPion->SetLineWidth(4);
1120
1121 TF1* MomentumFunctionProton = (TF1*) BetaGammaFunctionProton->Clone("MomentumFunctionProton");
1122 MomentumFunctionProton->SetParameter(2, m_ProtonPDGMass);
1123 MomentumFunctionProton->SetRange(0.01, 5.5);
1124 MomentumFunctionProton->SetLineColor(kRed);
1125 MomentumFunctionProton->SetLineWidth(4);
1126
1127 TF1* MomentumFunctionKaon = (TF1*) BetaGammaFunctionKaon->Clone("MomentumFunctionKaon");
1128 MomentumFunctionKaon->SetParameter(2, m_KaonPDGMass);
1129 MomentumFunctionKaon->SetRange(0.01, 5.5);
1130 MomentumFunctionKaon->SetLineColor(kRed);
1131 MomentumFunctionKaon->SetLineWidth(4);
1132
1133 gStyle->SetOptFit(1111);
1134 std::unique_ptr<TCanvas> CanvasOverlays(new TCanvas("CanvasOverlays", "overlays", 1300, 1000));
1135 CanvasOverlays->Divide(2, 2);
1136 CanvasOverlays->cd(1); Electron2DHistogram->Draw(); MomentumFunctionElectron->Draw("SAME");
1137 CanvasOverlays->cd(2); Pion2DHistogram->Draw(); MomentumFunctionPion->Draw("SAME");
1138 CanvasOverlays->cd(3); Kaon2DHistogram->Draw(); MomentumFunctionKaon->Draw("SAME");
1139 CanvasOverlays->cd(4); Proton2DHistogram->Draw(); MomentumFunctionProton->Draw("SAME");
1140 CanvasOverlays->Print("SVDdEdxOverlaysFitsHistos.pdf");
1141
1142 TF1* MomentumFunctionDeuteron = (TF1*) BetaGammaFunctionProton->Clone("MomentumFunctionDeuteron");
1143 MomentumFunctionDeuteron->SetParameter(2, m_DeuteronPDGMass);
1144 MomentumFunctionDeuteron->SetRange(0.01, 5.5);
1145 MomentumFunctionDeuteron->SetLineColor(kRed);
1146
1147 TF1* MomentumFunctionMuon = (TF1*) BetaGammaFunctionPion->Clone("MomentumFunctionMuon");
1148 MomentumFunctionMuon->SetParameter(2, m_MuonPDGMass);
1149 MomentumFunctionMuon->SetRange(0.01, 5.5);
1150 MomentumFunctionMuon->SetLineColor(kRed);
1151
1152// overlay all fits in one plot
1153
1154 std::unique_ptr<TCanvas> OverlayAllTracksCanvas(new TCanvas("OverlayAllTracksCanvas", "The Ultimate Plot", 10, 10, 1000, 700));
1155
1156 TH2F* AllTracksHistogram = new TH2F("AllTracksHistogram", "AllTracksHistogram;Momentum [GeV/c];dEdx [arb. units]", 1000, 0.05, 5,
1157 1000, 2.e5, 6.e6);
1158
1159 ttreeGeneric->Draw("TrackSVDdEdx:TrackSVDdEdxTrackMomentum>>AllTracksHistogram", "TracknSVDHits>7", "goff");
1160 AllTracksHistogram->Draw("COLZ");
1161 AllTracksHistogram->GetXaxis()->SetTitle("Momentum [GeV/c]");
1162 AllTracksHistogram->GetYaxis()->SetTitle("dE/dx [arbitrary units]");
1163 MomentumFunctionElectron->Draw("SAME");
1164 MomentumFunctionMuon->Draw("SAME");
1165 MomentumFunctionPion->Draw("SAME");
1166 MomentumFunctionKaon->Draw("SAME");
1167 MomentumFunctionProton->Draw("SAME");
1168 MomentumFunctionDeuteron->Draw("SAME");
1169 OverlayAllTracksCanvas->SetLogx();
1170 OverlayAllTracksCanvas->SetLogz();
1171
1172 OverlayAllTracksCanvas->Print("SVDdEdxAllTracksWithFits.pdf");
1173 TFile OverlayAllTracksPlotFile("SVDdEdxCalibrationOverlayAllTracks.root", "RECREATE");
1174 AllTracksHistogram->Write();
1175 MomentumFunctionElectron->Write();
1176 MomentumFunctionMuon->Write();
1177 MomentumFunctionPion->Write();
1178 MomentumFunctionKaon->Write();
1179 MomentumFunctionProton->Write();
1180 MomentumFunctionDeuteron->Write();
1181 OverlayAllTracksCanvas->Write();
1182 OverlayAllTracksPlotFile.Close();
1183
1184
1185// resolution studies //
1186
1187// 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.
1188 double PionRangeMin = 0.6;
1189 double PionRangeMax = 1.;
1190 double KaonRangeMin = 1.9;
1191 double KaonRangeMax = 3;
1192 double ElectronRangeMin = 1.;
1193 double ElectronRangeMax = 1.4;
1194
1195 auto PionResolutionHistogram = Pion2DHistogram->ProjectionY("PionResolutionHistogram",
1196 Pion2DHistogram->GetXaxis()->FindBin(PionRangeMin),
1197 Pion2DHistogram->GetXaxis()->FindBin(PionRangeMax));
1198 auto ElectronResolutionHistogram = Electron2DHistogram->ProjectionY("ElectronResolutionHistogram",
1199 Electron2DHistogram->GetXaxis()->FindBin(ElectronRangeMin), Electron2DHistogram->GetXaxis()->FindBin(ElectronRangeMax));
1200 auto KaonResolutionHistogram = Kaon2DHistogram->ProjectionY("KaonResolutionHistogram",
1201 Kaon2DHistogram->GetXaxis()->FindBin(KaonRangeMin),
1202 Kaon2DHistogram->GetXaxis()->FindBin(KaonRangeMax));
1203// for protons, there is not enough data in the flat range.
1204
1205
1206 TF1* PionResolutionFunction = new TF1("PionResolutionFunction",
1207 "[0]*TMath::Landau(x, [1], [1]*[2])*TMath::Gaus(x, [1], [1]*[2]*[4]) + [3]*TMath::Gaus(x, [1], [1]*[2]*[5])", 100e3, 1500e3);
1208// parameter [1] is the mean of the Landau
1209// parameter [2] is the relative resolution (w.r.t. the mean) of the Landau
1210// parameters [4]-[5] are the relative resolution of Gauss contributions w.r.t. that of Landau
1211// parameters [0] and [3] are fractions of the two components
1212 PionResolutionFunction->SetParameters(1, 6.e5, 0.1, 0.5, 2, 1);
1213 PionResolutionFunction->SetParLimits(0, 0, 500);
1214 PionResolutionFunction->SetParLimits(1, 3.e5, 8.e5);
1215 PionResolutionFunction->SetParLimits(2, 0, 1);
1216 PionResolutionFunction->SetParLimits(3, 0, 500);
1217 PionResolutionFunction->SetParLimits(4, 0, 7);
1218 PionResolutionFunction->SetParLimits(5, 1, 7);
1219 PionResolutionFunction->SetNpx(1000);
1220 auto FitResultResolutionPion = PionResolutionHistogram->Fit(PionResolutionFunction, "RSI");
1221
1222 B2INFO("relative resolution for pions: " << PionResolutionFunction->GetParameter(2));
1223 B2INFO("resolution for pions: fit status" << FitResultResolutionPion->Status());
1224
1225 TF1* KaonResolutionFunction = new TF1("KaonResolutionFunction",
1226 "[0]*TMath::Landau(x, [1], [1]*[2])*TMath::Gaus(x, [1], [1]*[2]*[4]) + [3]*TMath::Gaus(x, [1], [1]*[2]*[5])", 100e3, 1500e3);
1227
1228
1229 KaonResolutionFunction->SetParameters(1, 6.e5, 0.1, 0.5, 2, 1);
1230 KaonResolutionFunction->SetParLimits(0, 0, 500);
1231 KaonResolutionFunction->SetParLimits(1, 3.e5, 8.e5);
1232 KaonResolutionFunction->SetParLimits(2, 0, 1);
1233 KaonResolutionFunction->SetParLimits(3, 0, 500);
1234 KaonResolutionFunction->SetParLimits(4, 0, 7);
1235 KaonResolutionFunction->SetParLimits(5, 1, 7);
1236 KaonResolutionFunction->SetNpx(1000);
1237 auto FitResultResolutionKaon = KaonResolutionHistogram->Fit(KaonResolutionFunction, "RSI");
1238
1239 B2INFO("relative resolution for kaons: " << KaonResolutionFunction->GetParameter(2));
1240 B2INFO("resolution for kaons: fit status" << FitResultResolutionKaon->Status());
1241
1242 if ((FitResultResolutionKaon->Status() > 1)
1243 && (FitResultResolutionPion->Status() <= 1)) KaonResolutionFunction = (TF1*)
1244 PionResolutionFunction->Clone("KaonResolutionFunction");
1245
1246
1247
1248
1249 TF1* ElectronResolutionFunction = new TF1("ElectronResolutionFunction",
1250 "[0]*TMath::Landau(x, [1], [1]*[2])*TMath::Gaus(x, [1], [1]*[2]*[4]) + [3]*TMath::Gaus(x, [1], [1]*[2]*[5])", 50e3, 1500e3);
1251
1252
1253 ElectronResolutionFunction->SetParameters(1, 6.e5, 0.1, 0.5, 2, 1);
1254 ElectronResolutionFunction->SetParLimits(0, 0, 500);
1255 ElectronResolutionFunction->SetParLimits(1, 3.e5, 8.e5);
1256 ElectronResolutionFunction->SetParLimits(2, 0, 1);
1257 ElectronResolutionFunction->SetParLimits(3, 0, 500);
1258 ElectronResolutionFunction->SetParLimits(4, 0, 7);
1259 ElectronResolutionFunction->SetParLimits(5, 1, 7);
1260 ElectronResolutionFunction->SetNpx(1000);
1261 auto FitResultResolutionElectron = ElectronResolutionHistogram->Fit(ElectronResolutionFunction, "RSI");
1262
1263 B2INFO("relative resolution for electrons: " << ElectronResolutionFunction->GetParameter(2));
1264 B2INFO("resolution for electrons: fit status" << FitResultResolutionElectron->Status());
1265
1266 // plot all the resolution fits
1267 if (m_isMakePlots) {
1268 TCanvas* CanvasResolutions = new TCanvas("CanvasResolutions", "Resolutions", 1200, 650);
1269 CanvasResolutions->Divide(3, 1);
1270 CanvasResolutions->cd(1); PionResolutionHistogram->Draw();
1271 CanvasResolutions->cd(2); KaonResolutionHistogram->Draw();
1272 CanvasResolutions->cd(3); ElectronResolutionHistogram->Draw();
1273
1274 CanvasResolutions->Print("SVDdEdxResolutions.pdf");
1275 TFile OverlayResolutionsPlotFile("SVDdEdxCalibrationResolutions.root", "RECREATE");
1276 PionResolutionHistogram->Write();
1277 KaonResolutionHistogram->Write();
1278 ElectronResolutionHistogram->Write();
1279 CanvasResolutions->Write();
1280 OverlayResolutionsPlotFile.Close();
1281 }
1282
1283// evaluate the bias correction:
1284// difference between the MomentumFunction prediction and the mean of the resolution function in the flat part
1285// it should be of the order -1e4, i.e. about -1% of the absolute dEdx value
1286 double BiasCorrectionPion = PionResolutionFunction->GetParameter(1) - MomentumFunctionPion->Eval((
1287 PionRangeMax + PionRangeMin) / 2.);
1288 B2INFO("BiasCorrectionPion = " << BiasCorrectionPion);
1289
1290// generate a new pion payload using the MomentumFunctionPion, PionResolutionFunction and the bias correction
1291 TH2F* Pion2DHistogramNew = PrepareNewHistogram(Pion2DHistogram, Form("%sNew", Pion2DHistogram->GetName()), MomentumFunctionPion,
1292 PionResolutionFunction, BiasCorrectionPion);
1293
1294// sanity check: residual between the generated distribution and the data one
1295 TH2F* Pion2DHistogramResidual = (TH2F*) Pion2DHistogram->Clone("Pion2DHistogramResidual");
1296 Pion2DHistogramResidual->Add(Pion2DHistogramNew, Pion2DHistogram, 1, -1);
1297 Pion2DHistogramResidual->SetMinimum(-0.15);
1298 Pion2DHistogramResidual->SetMaximum(0.15);
1299
1300 // repeat, for kaons
1301 double BiasCorrectionKaon = KaonResolutionFunction->GetParameter(1) - MomentumFunctionKaon->Eval((
1302 KaonRangeMax + KaonRangeMin) / 2.);
1303 B2INFO("BiasCorrectionKaon = " << BiasCorrectionKaon);
1304
1305 // for protons, we compare the flat part of the MomentumFunction (~3 GeV) with the mean of the kaon resolution function
1306 // as there's not enough stats in the flat part to extract proton resolution from data
1307 double BiasCorrectionProton = KaonResolutionFunction->GetParameter(1) - MomentumFunctionProton->Eval(3.);
1308 B2INFO("BiasCorrectionProton = " << BiasCorrectionProton);
1309
1310 if ((BiasCorrectionProton / BiasCorrectionKaon) > 1.5) BiasCorrectionProton =
1311 BiasCorrectionKaon; // probably something went wrong due to low statistics
1312
1313 // back to kaons: generate a new payload
1314 TH2F* Kaon2DHistogramNew = PrepareNewHistogram(Kaon2DHistogram, Form("%sNew", Kaon2DHistogram->GetName()), MomentumFunctionKaon,
1315 KaonResolutionFunction, BiasCorrectionKaon);
1316// residual generated - data for kaons
1317 TH2F* Kaon2DHistogramResidual = (TH2F*) Kaon2DHistogram->Clone("Kaon2DHistogramResidual");
1318 Kaon2DHistogramResidual->Add(Kaon2DHistogramNew, Kaon2DHistogram, 1, -1);
1319 Kaon2DHistogramResidual->SetMinimum(-0.15);
1320 Kaon2DHistogramResidual->SetMaximum(0.15);
1321
1322 // same for protons (we use the kaon resolution function as explained above)
1323 TH2F* Proton2DHistogramNew = PrepareNewHistogram(Proton2DHistogram, Form("%sNew", Proton2DHistogram->GetName()),
1324 MomentumFunctionProton,
1325 KaonResolutionFunction, BiasCorrectionProton);
1326
1327// residual for protons
1328 TH2F* Proton2DHistogramResidual = (TH2F*) Proton2DHistogram->Clone("Proton2DHistogramResidual");
1329 Proton2DHistogramResidual->Add(Proton2DHistogramNew, Proton2DHistogram, 1, -1);
1330 Proton2DHistogramResidual->SetMinimum(-0.15);
1331 Proton2DHistogramResidual->SetMaximum(0.15);
1332
1333// deuterons: same as protons, but use the MomentumFunctionDeuteron
1334 TH2F* Deuteron2DHistogramNew = PrepareNewHistogram(Proton2DHistogram, "Deuteron2DHistogramNew", MomentumFunctionDeuteron,
1335 KaonResolutionFunction,
1336 BiasCorrectionKaon);
1337 Deuteron2DHistogramNew->SetTitle("hist_d1_1000010020_trunc");
1338
1339// muons: same as pions, but use the MomentumFunctionMuon
1340 TH2F* Muon2DHistogramNew = PrepareNewHistogram(Pion2DHistogram, "Muon2DHistogramNew", MomentumFunctionMuon, PionResolutionFunction,
1341 BiasCorrectionPion);
1342 Muon2DHistogramNew->SetTitle("hist_d1_13_trunc");
1343
1344// same for electrons
1345 double BiasCorrectionElectron = ElectronResolutionFunction->GetParameter(1) - MomentumFunctionElectron->Eval((
1346 ElectronRangeMax + ElectronRangeMin) / 2.);
1347 B2INFO("BiasCorrectionElectron = " << BiasCorrectionElectron);
1348 TH2F* Electron2DHistogramNew = PrepareNewHistogram(Electron2DHistogram, Form("%sNew", Electron2DHistogram->GetName()),
1349 MomentumFunctionElectron,
1350 ElectronResolutionFunction, BiasCorrectionElectron);
1351
1352 TH2F* Electron2DHistogramResidual = (TH2F*) Electron2DHistogram->Clone("Electron2DHistogramResidual");
1353 Electron2DHistogramResidual->Add(Electron2DHistogramNew, Electron2DHistogram, 1, -1);
1354 Electron2DHistogramResidual->SetMinimum(-0.15);
1355 Electron2DHistogramResidual->SetMaximum(0.15);
1356
1357 Electron2DHistogramNew->SetName("Electron2DHistogramNew");
1358 Muon2DHistogramNew->SetName("Muon2DHistogramNew");
1359 Pion2DHistogramNew->SetName("Pion2DHistogramNew");
1360 Kaon2DHistogramNew->SetName("Kaon2DHistogramNew");
1361 Proton2DHistogramNew->SetName("Proton2DHistogramNew");
1362 Deuteron2DHistogramNew->SetName("Deuteron2DHistogramNew");
1363
1364// plot the summary of all the distributions
1365 if (m_isMakePlots) {
1366 TCanvas* CanvasSummaryGenerated = new TCanvas("CanvasSummaryGenerated", "Generated payloads", 1700, 850);
1367 CanvasSummaryGenerated->Divide(3, 2);
1368 CanvasSummaryGenerated->cd(1); Electron2DHistogramNew->Draw("COLZ");
1369 CanvasSummaryGenerated->cd(2); Muon2DHistogramNew->Draw("COLZ");
1370 CanvasSummaryGenerated->cd(3); Pion2DHistogramNew->Draw("COLZ");
1371 CanvasSummaryGenerated->cd(4); Kaon2DHistogramNew->Draw("COLZ");
1372 CanvasSummaryGenerated->cd(5); Proton2DHistogramNew->Draw("COLZ");
1373 CanvasSummaryGenerated->cd(6); Deuteron2DHistogramNew->Draw("COLZ");
1374
1375 CanvasSummaryGenerated->Print("SVDdEdxGeneratedPayloads.pdf");
1376 TFile SummaryGeneratedPlotFile("SVDdEdxCalibrationSummaryGenerated.root", "RECREATE");
1377 Electron2DHistogramNew->Write();
1378 Muon2DHistogramNew->Write();
1379 Pion2DHistogramNew->Write();
1380 Kaon2DHistogramNew->Write();
1381 Proton2DHistogramNew->Write();
1382 Deuteron2DHistogramNew->Write();
1383 SummaryGeneratedPlotFile.Close();
1384
1385
1386 TCanvas* CanvasSummaryData = new TCanvas("CanvasSummaryData", "Data distributions", 1700, 850);
1387 CanvasSummaryData->Divide(3, 2);
1388 CanvasSummaryData->cd(1); Electron2DHistogram->Draw("COLZ");
1389 CanvasSummaryData->cd(3); Pion2DHistogram->Draw("COLZ");
1390 CanvasSummaryData->cd(4); Kaon2DHistogram->Draw("COLZ");
1391 CanvasSummaryData->cd(5); Proton2DHistogram->Draw("COLZ");
1392
1393 CanvasSummaryData->Print("SVDdEdxDataDistributions.pdf");
1394 TFile SummaryDataPlotFile("SVDdEdxCalibrationSummaryData.root", "RECREATE");
1395 Electron2DHistogram->Write();
1396 Pion2DHistogram->Write();
1397 Kaon2DHistogram->Write();
1398 Proton2DHistogram->Write();
1399 SummaryDataPlotFile.Close();
1400
1401
1402 TCanvas* CanvasSummaryResiduals = new TCanvas("CanvasSummaryResiduals", "Residuals", 1700, 850);
1403 CanvasSummaryResiduals->Divide(3, 2);
1404 CanvasSummaryResiduals->cd(1); Electron2DHistogramResidual->Draw("COLZ");
1405 CanvasSummaryResiduals->cd(3); Pion2DHistogramResidual->Draw("COLZ");
1406 CanvasSummaryResiduals->cd(4); Kaon2DHistogramResidual->Draw("COLZ");
1407 CanvasSummaryResiduals->cd(5); Proton2DHistogramResidual->Draw("COLZ");
1408
1409
1410 CanvasSummaryResiduals->Print("SVDdEdxResiduals.pdf");
1411 TFile SummaryResidualsPlotFile("SVDdEdxCalibrationSummaryResiduals.root", "RECREATE");
1412 Electron2DHistogramResidual->Write();
1413 Pion2DHistogramResidual->Write();
1414 Kaon2DHistogramResidual->Write();
1415 Proton2DHistogramResidual->Write();
1416 SummaryResidualsPlotFile.Close();
1417 }
1418
1419
1420 // return all the generated payloads
1421 std::unique_ptr<TList> histList(new TList);
1422 histList->Add(Electron2DHistogramNew);
1423 histList->Add(Muon2DHistogramNew);
1424 histList->Add(Pion2DHistogramNew);
1425 histList->Add(Kaon2DHistogramNew);
1426 histList->Add(Proton2DHistogramNew);
1427 histList->Add(Deuteron2DHistogramNew);
1428
1429 return histList;
1430
1431}
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::PrepareProfile
TH1F * PrepareProfile(TH2F *DataHistogram, TString NewName)
Reimplement the Profile histogram calculation for a 2D histogram.
Definition SVDdEdxCalibrationAlgorithm.h:232
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:523
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::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:311
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:667
Belle2::SVDdEdxCalibrationAlgorithm::DstarMassFit
TTree * DstarMassFit(std::shared_ptr< TTree > preselTree)
Mass fit for D*->Dpi.
Definition SVDdEdxCalibrationAlgorithm.cc:387
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:763
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