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