Belle II Software development
HadronCalibration.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 <cdc/calibration/CDCdEdx/HadronCalibration.h>
10#include <Math/MinimizerOptions.h>
11
12using namespace Belle2;
13
14HadronCalibration::HadronCalibration() {}
15
16
17void HadronCalibration::fitBGCurve(const std::vector< std::string >& particles, const std::string& filename,
18 const std::string& paramfile,
19 const std::string& suffix, const bool makeIterationSummary)
20{
21 ROOT::Math::MinimizerOptions::SetDefaultMinimizer("Minuit");
22 // read in a file that contains fit results for bg bins
23 double bgmin1 = 0.25, bgmax1 = 5.1; //using until 0 --> 4.5
24 double bgmin2 = 3.9, bgmax2 = 15.0; //using till 4.5 --> 10
25 double bgmin3 = 7.5, bgmax3 = 5000; //using above 10 --> 6000 *use this range only
26
27 const int npart = 5;
28 TFile* infile = new TFile(filename.data());
29
30 CDCDedxMeanPred gpar;
31 gpar.setParameters(paramfile);
32 CDCDedxWidgetCurve* gc = new CDCDedxWidgetCurve();
33
34 // multigraphs to hold the curve and residual results
35 TGraphErrors part_dedxvsbg[npart];
36 TMultiGraph* gr_dedxvsbg = new TMultiGraph(Form("gr_dedxvsbg_%s", suffix.data()), ";#beta#gamma;dE/dx");
37 TGraphErrors part_dedxvsp[npart];
38 TMultiGraph* gr_dedxvsp = new TMultiGraph(Form("gr_dedxvsp_%s", suffix.data()), ";momentum(GeV/c);dE/dx");
39
40 // --------------------------------------------------
41 // FILL BG CURVE VALUES
42 // --------------------------------------------------
43 TLegend tleg(0.4, 0.50, 0.65, 0.85);
44 tleg.SetBorderSize(0);
45
46 TLegend tlegPart(0.60, 0.60, 0.90, 0.90);
47 tlegPart.SetBorderSize(0);
48
49 for (int i = 0; i < int(particles.size()); ++i) {
50
51 std::string particle = particles[i];
52
53 double mass = m_prep.getParticleMass(particle);
54 if (mass == 0.0) B2FATAL("Mass of particle " << particle.data() << " is zero");
55
56 if (!infile->GetListOfKeys()->Contains(particle.data())) continue;
57
58 TTree* hadron = static_cast<TTree*>(infile->Get(particle.data()));
59 B2INFO("HadronCalibration: reading " << particle.data() << " in file " << filename.data());
60
61 double dedx, dedxerr, bg; // dE/dx without electron saturation correction
62
63 hadron->SetBranchAddress("dedx", &dedx);
64 hadron->SetBranchAddress("dedxerr", &dedxerr);
65 hadron->SetBranchAddress("bg_avg", &bg);
66
67 for (int j = 0; j < hadron->GetEntries(); ++j) {
68
69 hadron->GetEvent(j);
70
71 part_dedxvsbg[i].SetPoint(j, bg, dedx);
72 part_dedxvsbg[i].SetPointError(j, 0, dedxerr);
73
74 part_dedxvsp[i].SetPoint(j, bg * mass, dedx);
75 part_dedxvsp[i].SetPointError(j, 0, dedxerr);
76 }
77
78 part_dedxvsbg[i].SetName(particle.data());
79 if (i == 4) setGraphStyle(part_dedxvsbg[i], i + 2);
80 else setGraphStyle(part_dedxvsbg[i], i + 1);
81 gr_dedxvsbg->Add(&part_dedxvsbg[i]);
82
83 if (i == 4) setGraphStyle(part_dedxvsp[i], i + 2);
84 else setGraphStyle(part_dedxvsp[i], i + 1);
85 gr_dedxvsp->Add(&part_dedxvsp[i]);
86
87 tleg.AddEntry(&part_dedxvsbg[i], particles[i].data(), "p");
88 tlegPart.AddEntry(&part_dedxvsbg[i], particles[i].data(), "p");
89
90 }
91
92 TMultiGraph* grcopy1 = static_cast<TMultiGraph*>(gr_dedxvsbg->Clone(Form("datapoints_%s", suffix.data())));
93
94 //Region 1
95 // cppcheck-suppress constParameterPointer ; ROOT fixes this signature
96 TF1* fdedx1 = new TF1("fdedx1", [gc](double * x, double * par) {
97 std::vector<double> parVec(par, par + 8); // Create a vector from the parameters
98 return gc->meanCurve(x, parVec); // Call the member function
99 }, bgmin1, bgmax1, 8, "WidgetCurve");
100
101 // cppcheck-suppress constParameterPointer ; ROOT fixes this signature
102 TF1* fdedx1Copy = new TF1("fdedx1Copy", [gc](double * x, double * par) {
103 std::vector<double> parVec(par, par + 8); // Create a vector from the parameters
104 return gc->meanCurve(x, parVec); // Call the member function
105 }, bgmin1, bgmax1, 8, "WidgetCurve");
106
107 fdedx1->FixParameter(0, 1);
108 fdedx1Copy->FixParameter(0, 1);
109 for (int i = 1; i < 8; ++i) {
110 fdedx1->SetParameter(i, gpar.getCurvePars(i - 1));
111 fdedx1Copy->SetParameter(i, gpar.getCurvePars(i - 1));
112 }
113
114 //Region 2
115 // cppcheck-suppress constParameterPointer ; ROOT fixes this signature
116 TF1* fdedx2 = new TF1("fdedx2", [gc](double * x, double * par) {
117 std::vector<double> parVec(par, par + 5); // Create a vector from the parameters
118 return gc->meanCurve(x, parVec); // Call the member function
119 }, bgmin2, bgmax2, 5, "WidgetCurve");
120
121 // cppcheck-suppress constParameterPointer ; ROOT fixes this signature
122 TF1* fdedx2Copy = new TF1("fdedx2Copy", [gc](double * x, double * par) {
123 std::vector<double> parVec(par, par + 5); // Create a vector from the parameters
124 return gc->meanCurve(x, parVec); // Call the member function
125 }, bgmin2, bgmax2, 5, "WidgetCurve");
126
127 fdedx2->FixParameter(0, 2);
128 fdedx2Copy->FixParameter(0, 2);
129 for (int i = 1; i < 5; ++i) {
130 fdedx2->SetParameter(i, gpar.getCurvePars(6 + i));
131 fdedx2Copy->SetParameter(i, gpar.getCurvePars(6 + i));
132 }
133
134 //Region 3
135 // cppcheck-suppress constParameterPointer ; ROOT fixes this signature
136 TF1* fdedx3 = new TF1("fdedx3", [gc](double * x, double * par) {
137 std::vector<double> parVec(par, par + 5); // Create a vector from the parameters
138 return gc->meanCurve(x, parVec); // Call the member function
139 }, bgmin3, bgmax3, 5, "WidgetCurve");
140
141 // cppcheck-suppress constParameterPointer ; ROOT fixes this signature
142 TF1* fdedx3Copy = new TF1("fdedx3Copy", [gc](double * x, double * par) {
143 std::vector<double> parVec(par, par + 5); // Create a vector from the parameters
144 return gc->meanCurve(x, parVec); // Call the member function
145 }, bgmin3, bgmax3, 5, "WidgetCurve");
146
147 fdedx3->FixParameter(0, 3);
148 fdedx3Copy->FixParameter(0, 3);
149
150 for (int i = 1; i < 5; ++i) {
151 fdedx3->SetParameter(i, gpar.getCurvePars(10 + i));
152 fdedx3Copy->SetParameter(i, gpar.getCurvePars(10 + i));
153 }
154
155 // --------------------------------------------------
156 // FIT BG CURVE
157 // --------------------------------------------------
158
159 // //Fitting part1 1/beta^2 region and new constants
160 // //----------------------------------------------
161
162 int stat1 = gr_dedxvsbg->Fit("fdedx1", "", "", bgmin1, bgmax1);
163 for (int i = 0; i < 50; ++i) {
164 if (stat1 == 0) {
165 B2INFO("\n\tPART-1 FIT STATUS is OK: irr # " << i);
166 break;
167 }
168 stat1 = gr_dedxvsbg->Fit("fdedx1", "", "", bgmin1, bgmax1);
169 }
170
171 // if the fit was successful, write out the updated parameters
172 if (stat1 == 0) {
173 B2INFO("\t--> HadronCalibration: ATTENTIONS: PART-1 FIT OK..updating parameters");
174 for (int i = 1; i < 8; ++i) {
175 B2INFO("\t" << i << ") Old = " << gpar.getCurvePars(i - 1) << " --> New = " << fdedx1->GetParameter(i));
176 gpar.setCurvePars(i - 1, fdedx1->GetParameter(i));
177 }
178 } else B2INFO("\t--> HadronCalibration: WARNING: PART-1 FIT FAILED...");
179
180
181 //Fitting part2 min ionisation region and new constants
182 //----------------------------------------------
183 int stat2 = gr_dedxvsbg->Fit("fdedx2", "FQR", "", bgmin2, bgmax2);
184 for (int i = 0; i < 50; ++i) {
185 if (stat2 == 0) {
186 B2INFO("\n\tPART-2 FIT STATUS is OK: irr # " << i);
187 break;
188 }
189 stat2 = gr_dedxvsbg->Fit("fdedx2", "FQR", "", bgmin2, bgmax2);
190 }
191 if (stat2 == 0) {
192 B2INFO("\t--> HadronCalibration: ATTENTIONS: PART-2 FIT OK..updating parameters");
193 for (int i = 1; i < 5; ++i) {
194 B2INFO("\t" << i << ") Old = " << gpar.getCurvePars(6 + i) << " --> New = " << fdedx2->GetParameter(i));
195 gpar.setCurvePars(6 + i, fdedx2->GetParameter(i));
196 }
197 } else B2INFO("\t--> HadronCalibration: WARNING: PART-2 FIT FAILED...");
198
199
200 // //Fitting part3 relativistic region and new constants
201 // //----------------------------------------------
202 int stat3 = gr_dedxvsbg->Fit("fdedx3", "FQR", "", bgmin3, bgmax3);
203 for (int i = 0; i < 50; ++i) {
204 if (stat3 == 0) {
205 B2INFO("\n\tPART-3 FIT STATUS is OK: irr # " << i);
206 break;
207 }
208 stat3 = gr_dedxvsbg->Fit("fdedx3", "FQR", "", bgmin3, bgmax3);
209 }
210
211 if (stat3 == 0) {
212 B2INFO("\t--> HadronCalibration: ATTENTIONS: PART-3 FIT OK..updating parameters");
213 for (int i = 1; i < 5; ++i) {
214 B2INFO("\t" << i << ") Old = " << gpar.getCurvePars(10 + i) << " --> New = " << fdedx3->GetParameter(i));
215 gpar.setCurvePars(10 + i, fdedx3->GetParameter(i));
216 }
217 } else B2INFO("\t--> HadronCalibration: WARNING: PART-3 FIT FAILED...");
218
219 // //Plot without fitting (old fits + data points)
220 if (makeIterationSummary) {
221 TLine bgline1(4.5, 0.50, 4.5, 1.20);
222 bgline1.SetLineStyle(kDashed);
223 bgline1.SetLineColor(kGray);
224 TLine bgline2(10.0, 0.50, 10.0, 1.20);
225 bgline2.SetLineStyle(kDashed);
226 bgline2.SetLineColor(kGray);
227 TLine dedxline1(0.75, 1.0, 1000, 1.0);
228 dedxline1.SetLineStyle(kDashed);
229 dedxline1.SetLineColor(kGray);
230
231 setFitterStyle(fdedx1, 2, 1);
232 setFitterStyle(fdedx2, 5, 1);
233 setFitterStyle(fdedx3, 8, 1);
234
235 setFitterStyle(fdedx1Copy, 13, 6);
236 setFitterStyle(fdedx2Copy, 4, 6);
237 setFitterStyle(fdedx3Copy, 1, 6);
238
239 TCanvas bgcurvecan(Form("bgcurvecan_%s", suffix.data()), "bg curve and fitting", 1400, 800);
240 bgcurvecan.Divide(3, 2);
241 for (int i = 0; i < 6; i++) {
242 TMultiGraph* grcopy = static_cast<TMultiGraph*>(grcopy1->Clone(Form("datapoints_%s_%d", suffix.data(), i)));
243
244 bgcurvecan.cd(i + 1);
245 gPad->cd();
246 if (i == 0 || i == 3) gPad->SetLogy();
247 gPad->SetLogx();
248 grcopy->GetListOfGraphs()->SetDrawOption("AXIS");
249 grcopy->Draw("A*");
250 gPad->Modified(); gPad->Update();
251
252 if (i == 0 || i == 3) {
253 grcopy->GetXaxis()->SetLimits(0.10, 14500);
254 grcopy->SetMinimum(0.50);
255 grcopy->SetMaximum(50.0);
256 } else if (i == 1 || i == 4) {
257 grcopy->GetXaxis()->SetLimits(0.75, 100);
258 grcopy->SetMinimum(0.020);
259 grcopy->SetMaximum(3.0);
260 } else if (i == 2 || i == 5) {
261 grcopy->GetXaxis()->SetLimits(0.75, 50);
262 grcopy->SetMinimum(0.50);
263 grcopy->SetMaximum(1.20);
264 }
265 gPad->Modified(); gPad->Update();
266 TPaveText* ptold = new TPaveText(.35, .40, .75, .50, "blNDC");
267 ptold->SetBorderSize(0);
268 ptold->SetFillStyle(0);
269 if (i < 3) {
270 fdedx1Copy->Draw("same");
271 fdedx2Copy->Draw("same");
272 fdedx3Copy->Draw("same");
273 ptold->AddText("Old parameters");
274 } else {
275 fdedx1->Draw("same");
276 fdedx2->Draw("same");
277 fdedx3->Draw("same");
278 ptold->AddText("New parameters");
279
280 }
281
282 bgline1.Draw("same");
283 bgline2.Draw("same");
284 dedxline1.Draw("same");
285
286 if (i == 0) {
287 tleg.AddEntry(fdedx1Copy, "Pub-Fit: 1", "f1");
288 tleg.AddEntry(fdedx2Copy, "Pub-Fit: 2", "f1");
289 tleg.AddEntry(fdedx3Copy, "Pub-Fit: 3", "f1");
290 tleg.Draw("same");
291 }
292 if (i == 0 || i == 3) ptold->Draw("same");
293 }
294
295 bgcurvecan.SaveAs(Form("plots/HadronCal/BGfits/bgcurve_vsfits_%s.pdf", suffix.data()));
296
297 TCanvas bgcurveraw(Form("bgcurveraw_%s", suffix.data()), "bg curvs", 600, 600);
298 gPad->SetLogy();
299 gPad->SetLogx();
300 grcopy1->GetListOfGraphs()->SetDrawOption("AXIS");
301 grcopy1->Draw("A*");
302 gPad->Modified(); gPad->Update();
303 tlegPart.Draw("same");
304 fdedx1->Draw("same");
305 fdedx2->Draw("same");
306 fdedx3->Draw("same");
307 bgcurveraw.SaveAs(Form("plots/HadronCal/BGfits/bgcurve_raw_%s.root", suffix.data()));
308 bgcurveraw.SaveAs(Form("plots/HadronCal/BGfits/bgcurve_raw_%s.pdf", suffix.data()));
309
310 bgcurveraw.cd();
311 gPad->SetLogy();
312 gPad->SetLogx();
313 gr_dedxvsp->GetListOfGraphs()->SetDrawOption("AXIS");
314 gr_dedxvsp->Draw("A*");
315 gPad->Modified(); gPad->Update();
316 tlegPart.Draw("same");
317 bgcurveraw.SaveAs(Form("plots/HadronCal/BGfits/dedx_vs_mom_raw_%s.pdf", suffix.data()));
318 }
319
320 double func1a = fdedx1->Eval(4.5);
321 double func2a = fdedx2->Eval(4.5);
322 double func2b = fdedx2->Eval(10);
323 double func3a = fdedx3->Eval(10);
324 double diffval1 = 100 * abs(func1a - func2a) / func2a;
325 double diffval2 = 100 * abs(func2b - func3a) / func3a;
326 B2INFO("\t\n FIT Constraint for 1/beta^2 region (bg = 4.5): func1 --> " << func1a << ", func2 --> " << func2a <<
327 ", diff in % = " << diffval1);
328 B2INFO("\t\n FIT Constraint for 1/beta^2 region (bg = 10.): func1 --> " << func2b << ", func2 --> " << func3a <<
329 ", diff in % = " << diffval2);
330
331
332 // --------------------------------------------------
333 // GET RESIDUALS AND CHIS
334 // --------------------------------------------------
335 if (makeIterationSummary) {
336 TMultiGraph* fit_bgratio = new TMultiGraph(Form("fit_bgratio_%s", suffix.data()), ";#beta#gamma;ratio");
337 TGraph part_bgfit_ratio[npart];
338
339 TMultiGraph* fit_residual = new TMultiGraph(Form("fit_residual_%s", suffix.data()), ";#beta#gamma;residual");
340 TGraph part_bgfit_residual[npart];
341
342 double A = 4.5, B = 10;
343 int respoint = 1;
344 double rmin = 1.0, rmax = 1.0;
345
346 for (int i = 0; i < npart; ++i) {
347
348 for (int j = 0; j < part_dedxvsbg[i].GetN(); ++j) {
349
350 double x, y, fit;
351 part_dedxvsbg[i].GetPoint(j, x, y);
352
353 if (y == 0) continue;
354
355 if (x < A)
356 fit = fdedx1->Eval(x);
357 else if (x < B)
358 fit = fdedx2->Eval(x);
359 else
360 fit = fdedx3->Eval(x);
361
362 // the curve is just 1 for electrons...
363 if (npart == 4) fit = 1.0;
364 if (x < 2000) {
365 part_bgfit_ratio[i].SetPoint(respoint++, x, fit / y);
366 part_bgfit_residual[i].SetPoint(respoint++, x, fit - y);
367 }
368
369 if (fit / y < rmin) rmin = fit / y;
370 else if (fit / y > rmax) rmax = fit / y;
371
372 }
373
374 part_bgfit_ratio[i].SetMarkerSize(0.50);
375 part_bgfit_ratio[i].SetMarkerStyle(4);
376 part_bgfit_ratio[i].SetMarkerColor(i + 1);
377 if (i == 4) part_bgfit_ratio[i].SetMarkerColor(i + 2);
378 if (i <= 3)fit_bgratio->Add(&part_bgfit_ratio[i]);
379
380 part_bgfit_residual[i].SetMarkerSize(0.50);
381 part_bgfit_residual[i].SetMarkerStyle(4);
382 part_bgfit_residual[i].SetMarkerColor(i + 1);
383 if (i == 4)part_bgfit_residual[i].SetMarkerColor(i + 2);
384 if (i <= 3)fit_residual->Add(&part_bgfit_residual[i]);
385 }
386
387 fit_bgratio->SetMinimum(rmin * 0.97);
388 fit_bgratio->SetMaximum(rmax * 1.03);
389
390 TCanvas* bgfitratiocan = new TCanvas(Form("bgfitratiocan_%s", suffix.data()), "dE/dx fit residual", 450, 350);
391 bgfitratiocan->cd()->SetLogx();
392 bgfitratiocan->cd()->SetGridy();
393 fit_bgratio->Draw("AP");
394 tlegPart.Draw("same");
395 bgfitratiocan->SaveAs(Form("plots/HadronCal/BGfits/bgfit_ratios_%s.pdf", suffix.data()));
396 delete bgfitratiocan;
397
398 fit_residual->SetMinimum(-0.12);
399 fit_residual->SetMaximum(+0.12);
400
401 TCanvas* bgrescan = new TCanvas(Form("bgrescan_%s", suffix.data()), "dE/dx fit residual", 450, 350);
402 bgrescan->cd()->SetLogx();
403 bgrescan->cd()->SetGridy();
404 fit_residual->Draw("AP");
405 tlegPart.Draw("same");
406 bgrescan->SaveAs(Form("plots/HadronCal/BGfits/bgfit_residual_%s.pdf", suffix.data()));
407 delete bgrescan;
408 }
409 // write out the (possibly) updated parameters to file
410 gpar.printParameters("parameters.bgcurve.fit"); //Creating new file with new parameters
411
412 delete gr_dedxvsbg;
413 delete gr_dedxvsp;
414
415 delete fdedx1;
416 delete fdedx2;
417 delete fdedx3;
418 delete fdedx1Copy;
419 delete fdedx2Copy;
420 delete fdedx3Copy;
421
422 infile->Close();
423}
424
425void HadronCalibration::plotBGMonitoring(const std::vector< std::string >& particles, const std::string& filename,
426 const std::string& suffix)
427{
428 //1. chi-mean vs bg
429 std::string title = "chi-fit-means of different particles;#beta#gamma;#chi(#mu)";
430 std::string sname = Form("gr_chimean_vs_bg_%s", suffix.data());
431 plotMonitoring(particles, filename, sname, title, "bg", "chi");
432
433 //2. Sigma vs bg
434 title = "chi-fit-width of different particles;#beta#gamma;#chi(#sigma)";
435 sname = Form("gr_chiwidth_vs_bg_%s", suffix.data());
436 plotMonitoring(particles, filename, sname, title, "bg", "sigma");
437
438 //3. chi mean vs p
439 title = "chi-fit-means of different particles: latest curve+sigma pars;p(GeV/c);#chi(#mu)";
440 sname = Form("gr_chimean_vs_mom_%s", suffix.data());
441 plotMonitoring(particles, filename, sname, title, "mom", "chi");
442
443 //4. Sigma vs p
444 title = "chi-fit-width of different particles: w/ latest curve+sigma pars;p(GeV/c);#chi(#sigma)";
445 sname = Form("gr_chiwidth_vs_mom_%s", suffix.data());
446 plotMonitoring(particles, filename, sname, title, "mom", "sigma");
447
448 std::string svar = "ler";
449 for (int ir = 0; ir < 2; ir++) {
450 if (ir == 1) svar = "her";
451
452 //5. chi mean vs injection time
453 title = Form("#chi mean vs injection time, %s;injection time;#chi_{#mu}", svar.data());
454 sname = Form("gr_chimean_vs_inj_%s_%s", suffix.data(), svar.data());
455 plotMonitoring(particles, filename, sname, title, svar, "chi");
456
457 //5. chi sigma vs injection time
458 title = Form("#chi sigma vs injection time, %s;injection time;#chi_{#sigma}", svar.data());
459 sname = Form("gr_chiwidth_vs_inj_%s_%s", suffix.data(), svar.data());
460 plotMonitoring(particles, filename, sname, title, svar, "sigma");
461 }
462}
463
464void HadronCalibration::plotMonitoring(const std::vector< std::string >& particles, const std::string& filename,
465 const std::string& sname,
466 const std::string& title,
467 const std::string& sxvar, const std::string& syvar)
468{
469
470 const int npart = int(particles.size());
471
472 // read in a file that contains fit results for bg bins
473 TFile* infile = new TFile(filename.data());
474
475 // multigraphs to hold the curve and residual results
476 // TGraphErrors grchim[npart];
477 std::vector<TGraphErrors> grchim(npart);
478
479 TMultiGraph* gr_var = new TMultiGraph(Form("%s", sname.data()), "");
480
481 TLegend tlegPart(0.75, 0.65, 0.90, 0.90);
482 tlegPart.SetBorderSize(0);
483
484 for (int i = 0; i < npart; ++i) {
485 std::string particle = particles[i];
486 double mass = m_prep.getParticleMass(particle);
487 if (mass == 0.0) B2FATAL("Mass of particle " << particle.data() << " is zero");
488
489 if (!infile->GetListOfKeys()->Contains(particle.data())) continue;
490 TTree* hadron;
491 if (sxvar == "bg" || sxvar == "mom") hadron = static_cast<TTree*>(infile->Get(particle.data()));
492 else hadron = static_cast<TTree*>(infile->Get(Form("%s_%s", particle.data(), sxvar.data())));
493
494 B2INFO("HadronCalibration: reading " << particle.data() << " in file " << filename.data());
495
496 double chimean, chisigma, avg, chimean_err, chisigmaerr;
497
498 hadron->SetBranchAddress("chimean", &chimean);
499 hadron->SetBranchAddress("chimean_err", &chimean_err);
500 hadron->SetBranchAddress("chisigma", &chisigma);
501 hadron->SetBranchAddress("chisigma_err", &chisigmaerr);
502 if (sxvar == "bg" || sxvar == "mom") hadron->SetBranchAddress("bg_avg", &avg);
503 else hadron->SetBranchAddress("inj_avg", &avg);
504
505 for (int j = 0; j < hadron->GetEntries(); ++j) {
506
507 hadron->GetEvent(j);
508 double var, xvar; //varerr
509 if (syvar == "chi") {var = chimean;} // varerr = chimean_err;}
510 else {var = chisigma;} // varerr = chisigmaerr;}
511 if (sxvar == "mom") xvar = avg * mass;
512 else xvar = avg;
513
514 grchim[i].SetPoint(j, xvar, var);
515 //grchim[i].SetPointError(j, 0, varerr);
516 }
517
518 tlegPart.AddEntry(&grchim[i], particles[i].data(), "p");
519
520 setGraphStyle(grchim[i], i + 1);
521 if (i == 4) grchim[i].SetMarkerColor(i + 2);
522 gr_var->Add(&grchim[i]);
523
524 }
525
526 TCanvas ctemp(Form("%s", sname.data()), Form("%s", sname.data()), 450, 350);
527 gPad->SetGridy();
528 if (sxvar == "bg" || sxvar == "mom") gPad->SetLogx();
529
530 if (syvar == "chi") { gr_var->SetMinimum(-1.0); gr_var->SetMaximum(+1.0);}
531 else { gr_var->SetMinimum(0.6); gr_var->SetMaximum(1.4); }
532
533 gr_var->SetTitle(Form("%s", title.data()));
534 gr_var->Draw("AP");
535 tlegPart.Draw("same");
536 ctemp.SaveAs(Form("plots/HadronCal/Monitoring/%s.pdf", sname.data()));
537 if (sxvar == "bg" || sxvar == "mom") {
538 gPad->SetLogx();
539 if (sxvar == "bg") gr_var->GetXaxis()->SetLimits(0.1, 20);
540 else gr_var->GetXaxis()->SetLimits(0.1, 1.0);
541 gr_var->Draw("AP");
542 tlegPart.Draw("same");
543 ctemp.SaveAs(Form("plots/HadronCal/Monitoring/%s_zoomed.pdf", sname.data()));
544 }
545 delete gr_var;
546}
547
548void HadronCalibration::fitSigmavsIonz(const std::vector< std::string >& particles, const std::string& filename,
549 const std::string& paramfile,
550 const std::string& suffix, const bool makeIterationSummary)
551{
552
553 // read in a file that contains fit results for bg bins
554 const int npart = int(particles.size());
555
556 CDCDedxSigmaPred gpar;
557 gpar.setParameters(paramfile);
558
559 TFile* infile = new TFile(filename.data());
560 std::vector<TGraphErrors> part_resovsdedx(npart);
561
562 TMultiGraph* gr_resovsdedx = new TMultiGraph("gr_resovsdedx", ";dedx;#sigma(ionz)");
563
564 TLegend tlegPart(0.72, 0.15, 0.88, 0.40);
565
566 for (int i = 0; i < int(particles.size()); ++i) {
567
568 std::string particle = particles[i];
569
570 double mass = m_prep.getParticleMass(particle);
571 if (mass == 0.0) B2FATAL("Mass of particle " << particle.data() << " is zero");
572
573 if (particle == "electron") continue;
574
575 if (!infile->GetListOfKeys()->Contains(particle.data())) continue;
576 TTree* hadron = static_cast<TTree*>(infile->Get(particle.data()));
577 B2INFO("\tHadronCalibration: reading " << particle << " in file " << filename.data());
578
579 double dedx, ionzres;
580
581 hadron->SetBranchAddress("dedx", &dedx);
582 hadron->SetBranchAddress("ionzres", &ionzres);
583
584 part_resovsdedx[i].SetName(particle.data());
585
586 for (int j = 0; j < hadron->GetEntries(); ++j) {
587 hadron->GetEvent(j);
588 part_resovsdedx[i].SetPoint(j, dedx, ionzres);
589 }
590
591 if (i == 4) setGraphStyle(part_resovsdedx[i], i + 2);
592 else setGraphStyle(part_resovsdedx[i], i + 1);
593 gr_resovsdedx->Add(&part_resovsdedx[i]);
594 tlegPart.AddEntry(&part_resovsdedx[i], particles[i].data(), "p");
595 }
596
597 gStyle->SetOptStat(0);
598 gStyle->SetStatY(0.9); // Set y-position (fraction of pad size)
599 gStyle->SetStatX(0.4); // Set x-position (fraction of pad size)
600 gStyle->SetStatW(0.15); // Set width of stat-box (fraction of pad size)
601 gStyle->SetStatH(0.15); // Set height of stat-box (fraction of pad size)
602
603 TF1* sigvsdedx = new TF1("sigvsdedx", "[0]+[1]*x", 0.50, 15.0);
604 sigvsdedx->SetParameter(0, gpar.getDedxPars(0));
605 sigvsdedx->SetParameter(1, gpar.getDedxPars(1));
606
607 TF1* sigvsdedxCopy = static_cast<TF1*>(sigvsdedx->Clone("sigvsdedxcopy"));
608 setFitterStyle(sigvsdedxCopy, 13, 6);
609
610 TCanvas sigcan("sigcan", " Reso(ionz) vs dE/dx", 820, 750);
611 sigcan.cd();
612 gr_resovsdedx->Draw("APE");
613 sigvsdedxCopy->Draw("same");
614 tlegPart.Draw("same");
615
616 int status = gr_resovsdedx->Fit("sigvsdedx", "MF", "", 0.50, 7.0);
617 for (int i = 0; i < 10; ++i) {
618 if (status == 0) break;
619 status = gr_resovsdedx->Fit("sigvsdedx", "", "", 0.50, 7.0);
620 }
621 if (makeIterationSummary) {
622 gr_resovsdedx->SetTitle(Form("%s (slope = %0.03f, const = %0.03f)", gr_resovsdedx->GetTitle(), sigvsdedx->GetParameter(1),
623 sigvsdedx->GetParameter(0)));
624 gPad->Modified(); gPad->Update();
625 sigcan.SaveAs(Form("plots/HadronCal/Resofits/sigma_vsionz_%s.pdf", suffix.data()));
626
627 gr_resovsdedx->GetXaxis()->SetLimits(0.2, 2.00);
628 gr_resovsdedx->GetHistogram()->SetMaximum(0.50);
629 gr_resovsdedx->GetHistogram()->SetMinimum(0.00);
630 gr_resovsdedx->Draw("APE");
631 sigvsdedxCopy->Draw("same");
632 tlegPart.Draw("same");
633 sigcan.SaveAs(Form("plots/HadronCal/Resofits/sigma_vsionz_zoomed_%s.pdf", suffix.data()));
634 gStyle->SetOptStat(11);
635 }
636 // if the fit was successful, save the updated parameters
637 if (status == 0) {
638 B2INFO("\tHadronCalibration: SigmavsdEdx FITs Ok. updating parameters");
639 for (int i = 0; i < 2; ++i) {
640 B2INFO("\t" << i << ") Old = " << gpar.getDedxPars(i) << " --> New = " << sigvsdedx->GetParameter(i));
641 gpar.setDedxPars(i, sigvsdedx->GetParameter(i));
642 }
643 } else B2INFO("\tHadronCalibration: WARNING: SigmavsdEdx FIT FAILED... \n \tHadronCalibration: Skipping parameters");
644
645 // write out the (possibly) updated parameters to file
646 gpar.printParameters("parameters.ionz.fit"); //Creating new file with new parameters
647 infile->Close();
648
649 delete gr_resovsdedx;
650 delete sigvsdedx;
651}
652
653void HadronCalibration::fitSigmaVsNHit(const std::vector< std::string >& particles, const std::string& filename,
654 const std::string& paramsigma,
655 const std::string& suffix, const bool makeIterationSummary)
656{
657
658 const double lowernhit = 7, uppernhit = 39;
659
660 CDCDedxSigmaPred sgpar;
661 sgpar.setParameters(paramsigma);
662
663 CDCDedxWidgetSigma* gs = new CDCDedxWidgetSigma();
664
665 // cppcheck-suppress constParameterPointer ; ROOT fixes this signature
666 TF1* fsigma = new TF1("fsigma", [gs](double * x, double * par) {
667 std::vector<double> parVec(par, par + 6); // Create a vector from the parameters
668 return gs->sigmaCurve(x, parVec); // Call the member function
669 }, lowernhit, uppernhit, 6, "CDCDedxWidgetSigma");
670
671 fsigma->FixParameter(0, 2);
672 for (int i = 1; i < 6; ++i) fsigma->SetParameter(i, sgpar.getNHitPars(i - 1));
673
674 TF1* fsigmacopy = static_cast<TF1*>(fsigma->Clone("fsigmacopy"));
675 setFitterStyle(fsigmacopy, 13, 6);
676
677 TFile* infile = new TFile(filename.data());
678
679 for (int ip = 0; ip < int(particles.size()); ++ip) {
680
681 std::string particle = particles[ip];
682 TGraphErrors gr_dedxvsbg;
683
684 if (!infile->GetListOfKeys()->Contains(Form("%s_nhit", particle.data()))) continue;
685 TTree* hadron = static_cast<TTree*>(infile->Get(Form("%s_nhit", particle.data())));
686 B2INFO("\tHadronCalibration: reading " << particle << " in file " << filename.data());
687
688 double avg, sigma, sigmaerr;
689
690 hadron->SetBranchAddress("avg", &avg);
691 hadron->SetBranchAddress("chisigma", &sigma);
692 hadron->SetBranchAddress("chisigma_err", &sigmaerr);
693
694 for (int j = 0; j < hadron->GetEntries(); ++j) {
695 hadron->GetEvent(j);
696 gr_dedxvsbg.SetPoint(j, avg, sigma);
697 gr_dedxvsbg.SetPointError(j, 0, sigmaerr);
698 }
699
700 // --------------------------------------------------
701 // FIT SIGMA VS NHIT CURVE
702 // --------------------------------------------------
703 gStyle->SetOptFit(0);
704 TCanvas sigvsnhitcan(Form("sigvsnhitcan_%s", suffix.data()), "#sigma vs. nHit", 600, 600);
705
706 gr_dedxvsbg.SetMarkerStyle(8);
707 gr_dedxvsbg.SetMarkerSize(0.5);
708 gr_dedxvsbg.SetMaximum(2.2);
709 gr_dedxvsbg.SetMinimum(0.0);
710 gr_dedxvsbg.SetTitle(Form("width of (dedx-pred)/(#sigma_{Cos * Ion * InjReso}) vs lNHitsUsed, %s ;lNHitsUsed; #sigma",
711 particle.data()));
712 gr_dedxvsbg.Draw("AP");
713 fsigmacopy->Draw("same");
714
715 TLegend tleg(0.4, 0.70, 0.65, 0.85);
716 tleg.SetBorderSize(0);
717 tleg.AddEntry(&gr_dedxvsbg, "Data points", "P");
718 tleg.AddEntry(fsigmacopy, "Public Fit", "f1");
719
720 if (particle == "muon") {
721
722 // if the fit succeeds, write out the new parameters
723 int status = gr_dedxvsbg.Fit("fsigma", "FR", "", lowernhit, uppernhit);
724
725 for (int i = 0; i < 20; ++i) {
726 if (status == 0) break;
727 status = gr_dedxvsbg.Fit("fsigma", "FR", "", lowernhit, uppernhit);
728 }
729
730 if (status == 0) {
731 B2INFO("\tHadronCalibration::fitSigmaVsNHit --> FIT OK..updating parameters");
732 for (int j = 1; j < 6; ++j) {
733 B2INFO("\t" << j << ") Old = " << sgpar.getNHitPars(j - 1) << " --> New = " << fsigma->GetParameter(j));
734 sgpar.setNHitPars(j - 1, fsigma->GetParameter(j));
735 }
736 } else B2INFO("\tHadronCalibration::fitSigmaVsNHit --> WARNING: FIT FAILED..: status = " << status);
737
738 // write out the (possibly) updated parameters to file
739 sgpar.printParameters("parameters.sigmanhit.fit");
740
741 tleg.AddEntry(fsigma, "New Fit", "f1");
742 }
743
744 else {
745 sgpar.setParameters("parameters.sigmanhit.fit");
746 for (int i = 1; i < 6; ++i) fsigma->SetParameter(i, sgpar.getNHitPars(i - 1));
747 fsigma->SetRange(7, 39);
748 fsigma->Draw("same");
749 tleg.AddEntry(fsigma, "New (param. from muon fit)", "f1");
750 }
751
752 tleg.Draw("same");
753
754 if (makeIterationSummary) sigvsnhitcan.SaveAs(Form("plots/HadronCal/Resofits/sigma_vsnhits_%s_%s.pdf", suffix.data(),
755 particle.data()));
756 }
757}
758
759void HadronCalibration::fitSigmaVsCos(const std::vector< std::string >& particles, const std::string& filename,
760 const std::string& paramsigma,
761 const std::string& suffix, const bool makeIterationSummary)
762{
763 double lowercos = -0.84, uppercos = 0.96;
764
765 CDCDedxSigmaPred gpar;
766 gpar.setParameters(paramsigma);
767
768 CDCDedxWidgetSigma* gs = new CDCDedxWidgetSigma();
769
770 // cppcheck-suppress constParameterPointer ; ROOT fixes this signature
771 TF1* total = new TF1("total", [gs](double * x, double * par) {
772 std::vector<double> parVec(par, par + 11); // Create a vector from the parameters
773 return gs->sigmaCurve(x, parVec); // Call the member function
774 }, lowercos, uppercos, 11, "CDCDedxWidgetSigma"); // 6 parameters for this example
775
776 for (int i = 0; i < 10; ++i) total->SetParameter(i + 1, gpar.getCosPars(i));
777 total->FixParameter(0, 3);
778 total->FixParameter(2, 0.0);
779
780 TF1* fsigmacopy = static_cast<TF1*>(total->Clone("fsigmacopy"));
781 setFitterStyle(fsigmacopy, 13, 6);
782
783 TFile* infile = new TFile(filename.data());
784
785 for (int ip = 0; ip < int(particles.size()); ++ip) {
786
787 std::string particle = particles[ip];
788 TGraphErrors gr_dedx;
789
790 if (!infile->GetListOfKeys()->Contains(Form("%s_costh", particle.data()))) continue;
791 TTree* hadron = static_cast<TTree*>(infile->Get(Form("%s_costh", particle.data())));
792 B2INFO("\tHadronCalibration: reading " << particle << " in file " << filename.data());
793
794 double avg, sigma, sigmaerr;
795
796 hadron->SetBranchAddress("avg", &avg);
797 hadron->SetBranchAddress("chisigma", &sigma);
798 hadron->SetBranchAddress("chisigma_err", &sigmaerr);
799
800 for (int j = 0; j < hadron->GetEntries(); ++j) {
801 hadron->GetEvent(j);
802 gr_dedx.SetPoint(j, avg, sigma);
803 gr_dedx.SetPointError(j, 0, sigmaerr);
804 }
805
806 // --------------------------------------------------
807 // FIT SIGMA VS COS CURVE
808 // --------------------------------------------------
809 gStyle->SetOptFit(0);
810
811 TCanvas sigvscos("sigvscos", "#sigma vs. cos(#theta)", 400, 400);
812
813 gr_dedx.SetMaximum(1.8);
814 gr_dedx.SetMinimum(0.4);
815 gr_dedx.SetTitle(Form("(dedx-pred)/(#sigma_{Nhit * Ion * InjReso}) vs cos(#theta), %s;cos(#theta);#sigma", particle.data()));
816 gr_dedx.Draw("AP");
817
818 fsigmacopy->Draw("same");
819
820 TLegend tleg(0.75, 0.75, 0.89, 0.89);
821 tleg.SetBorderSize(0);
822 tleg.AddEntry(&gr_dedx, "Data points", "P");
823 tleg.AddEntry(fsigmacopy, "Public Fit", "f1");
824
825 if (particle == "muon") {
826
827 int status = gr_dedx.Fit("total", "FR", "", lowercos, uppercos);
828
829 for (int i = 0; i < 20; ++i) {
830 if (status == 0) break;
831 status = gr_dedx.Fit("total", "FR", "", lowercos, uppercos);
832 }
833 if (status == 0) {
834 B2INFO("\tHadronCalibration::fitSigmaVsCos --> FIT OK..updating parameters");
835 for (int j = 1; j < 11; ++j) {
836 B2INFO("\t" << j - 1 << ") Old = " << gpar.getCosPars(j - 1) << " --> New = " << total->GetParameter(j));
837 gpar.setCosPars(j - 1, total->GetParameter(j));
838 }
839 } else B2INFO("\tHadronCalibration::fitSigmaVsCos --> WARNING: FIT FAILED..: status = " << status);
840
841 gpar.printParameters("parameters.sigmacos.fit");
842 tleg.AddEntry(total, "New Fit", "f1");
843
844 } else {
845 gpar.setParameters("parameters.sigmacos.fit");
846 for (int i = 1; i < 11; ++i) total->SetParameter(i, gpar.getCosPars(i - 1));
847 total->Draw("same");
848 tleg.AddEntry(total, "New (param. from muon fit)", "f1");
849
850 }
851
852 tleg.Draw("same");
853 if (makeIterationSummary) sigvscos.SaveAs(Form("plots/HadronCal/Resofits/sigma_vscos_%s_%s.pdf", suffix.data(), particle.data()));
854 }
855
856}
Belle2::CDCDedxMeanPred
Class to hold the prediction of mean as a function of beta-gamma (bg)
Definition CDCDedxMeanPred.h:31
Belle2::CDCDedxMeanPred::setCurvePars
void setCurvePars(int i, double val)
set the curve parameters
Definition CDCDedxMeanPred.h:63
Belle2::CDCDedxMeanPred::setParameters
void setParameters(const std::string &infile)
set the parameters from file
Definition CDCDedxMeanPred.cc:28
Belle2::CDCDedxMeanPred::getCurvePars
double getCurvePars(int i)
get the curve parameters
Definition CDCDedxMeanPred.h:58
Belle2::CDCDedxMeanPred::printParameters
void printParameters(const std::string &infile)
write the parameters in file
Definition CDCDedxMeanPred.cc:49
Belle2::CDCDedxSigmaPred
Class to hold the prediction of resolution depending dE/dx, nhit, and cos(theta)
Definition CDCDedxSigmaPred.h:32
Belle2::CDCDedxSigmaPred::setNHitPars
void setNHitPars(int i, double val)
set the nhit parameters
Definition CDCDedxSigmaPred.h:89
Belle2::CDCDedxSigmaPred::setParameters
void setParameters(const std::string &infile)
set the parameters from file
Definition CDCDedxSigmaPred.cc:27
Belle2::CDCDedxSigmaPred::getCosPars
double getCosPars(int i)
get the cos(theta) parameters
Definition CDCDedxSigmaPred.h:94
Belle2::CDCDedxSigmaPred::getDedxPars
double getDedxPars(int i)
get the dedx parameters
Definition CDCDedxSigmaPred.h:74
Belle2::CDCDedxSigmaPred::setDedxPars
void setDedxPars(int i, double val)
set the dedx parameters
Definition CDCDedxSigmaPred.h:79
Belle2::CDCDedxSigmaPred::setCosPars
void setCosPars(int i, double val)
set the cos(theta) parameters
Definition CDCDedxSigmaPred.h:99
Belle2::CDCDedxSigmaPred::getNHitPars
double getNHitPars(int i)
get the nhit parameters
Definition CDCDedxSigmaPred.h:84
Belle2::CDCDedxSigmaPred::printParameters
void printParameters(const std::string &infile)
write the parameters in file
Definition CDCDedxSigmaPred.cc:60
Belle2::CDCDedxWidgetCurve
Class to hold the beta-gamma (bg) mean function.
Definition CDCDedxWidgetCurve.h:23
Belle2::CDCDedxWidgetCurve::meanCurve
static double meanCurve(const double *x, const std::vector< double > &par)
calculate the predicted mean value as a function of beta-gamma (bg) this is done with a different fun...
Definition CDCDedxWidgetCurve.h:42
Belle2::CDCDedxWidgetSigma
Class to hold the beta-gamma (bg) resolution function.
Definition CDCDedxWidgetSigma.h:22
Belle2::CDCDedxWidgetSigma::sigmaCurve
static double sigmaCurve(const double *x, const std::vector< double > &par)
calculate the predicted sigma value as a function of beta-gamma (bg) this is done with a different fu...
Definition CDCDedxWidgetSigma.h:40
Belle2::HadronCalibration::plotBGMonitoring
void plotBGMonitoring(const std::vector< std::string > &particles, const std::string &filename, const std::string &suffix)
plots mean and width after fitting
Definition HadronCalibration.cc:425
Belle2::HadronCalibration::fitSigmaVsNHit
static void fitSigmaVsNHit(const std::vector< std::string > &particles, const std::string &filename, const std::string &paramsigma, const std::string &suffx, const bool makeIterationSummary)
fit sigma vs.
Definition HadronCalibration.cc:653
Belle2::HadronCalibration::plotMonitoring
void plotMonitoring(const std::vector< std::string > &particles, const std::string &filename, const std::string &sname, const std::string &title, const std::string &sx, const std::string &sy)
plots chi and width after fitting - main function
Definition HadronCalibration.cc:464
Belle2::HadronCalibration::fitSigmavsIonz
void fitSigmavsIonz(const std::vector< std::string > &particles, const std::string &filename, const std::string &paramfile, const std::string &suffix, const bool makeIterationSummary)
fit sigma vs.
Definition HadronCalibration.cc:548
Belle2::HadronCalibration::fitBGCurve
void fitBGCurve(const std::vector< std::string > &particles, const std::string &filename, const std::string &paramfile, const std::string &suffx, const bool makeIterationSummary)
fit the beta-gamma curve
Definition HadronCalibration.cc:17
Belle2::HadronCalibration::setFitterStyle
static void setFitterStyle(TF1 *&fitt, const int ic, const int il)
function to set fitter cosmetics
Definition HadronCalibration.h:109
Belle2::HadronCalibration::fitSigmaVsCos
static void fitSigmaVsCos(const std::vector< std::string > &particles, const std::string &filename, const std::string &paramfile, const std::string &suffx, const bool makeIterationSummary)
fit sigma vs.
Definition HadronCalibration.cc:759
Belle2::HadronCalibration::m_prep
HadronBgPrep m_prep
object for dE/dx to prepare sample
Definition HadronCalibration.h:118
Belle2::HadronCalibration::setGraphStyle
static void setGraphStyle(TGraphErrors &gr, const int ic)
function to set graph cosmetics
Definition HadronCalibration.h:99
Belle2::HadronCalibration::HadronCalibration
HadronCalibration()
Constructor: Sets the description, the properties and the parameters of the algorithm.
Definition HadronCalibration.cc:14
Belle2
Abstract base class for different kinds of events.
Definition MillepedeAlgorithm.h:17