9#include <cdc/calibration/CDCdEdx/CDCDedxInjectTimeAlgorithm.h>
11#include <cdc/utilities/CDCDedxMeanPred.h>
12#include <cdc/utilities/CDCDedxSigmaPred.h>
14#include <framework/gearbox/Const.h>
43 setDescription(
"A calibration algorithm for CDC dE/dx injection time gain and reso");
56 B2FATAL(
"There is no valid payload for Injection time");
62 double dedx = 0.0, injtime = 0.0, injring = 1.0, costh, mom;
64 ttree->SetBranchAddress(
"dedx", &dedx);
65 ttree->SetBranchAddress(
"injtime", &injtime);
66 ttree->SetBranchAddress(
"injring", &injring);
67 ttree->SetBranchAddress(
"costh", &costh);
68 ttree->SetBranchAddress(
"p", &mom);
69 ttree->SetBranchAddress(
"nhits", &nhits);
82 std::array<std::vector<TH1D*>, numdedx::nrings> hdedx, htime, hchi, hdedx_corr;
88 const double tzedges[4] = {0, 2.0e3, 0.1e6, 20e6};
89 std::array<std::array<TH1D*, 3>, numdedx::nrings> hztime;
93 for (
int i = 0; i < ttree->GetEntries(); ++i) {
96 if (dedx <= 0 || injtime < 0 || injring < 0)
continue;
103 if (injring > 0.5) wr = 1;
106 unsigned int tb = htimes->GetXaxis()->FindBin(injtime);
110 htimes->Fill(injtime);
111 if (injtime < tzedges[1]) hztime[wr][0]->Fill(injtime);
112 else if (injtime < tzedges[2]) hztime[wr][1]->Fill(injtime);
113 else hztime[wr][2]->Fill(injtime);
115 hdedx[wr][tb]->Fill(dedx);
116 htime[wr][tb]->Fill(injtime);
133 std::map<int, std::vector<double>> vmeans, vresos, vtimes;
136 for (
unsigned int ir = 0; ir <
c_rings; ir++) {
137 for (
unsigned int it = 0; it <
m_tbins; it++) {
138 double avgtime = htime[ir][it]->GetMean();
139 double avgtimeerr = htime[ir][it]->GetMeanError();
140 vtimes[ir * 2].push_back(avgtime);
141 vtimes[ir * 2 + 1].push_back(avgtimeerr);
149 std::map<int, std::vector<double>> vmeanscorr;
153 std::array<double, numdedx::nrings> scale;
154 std::map<int, std::vector<double>> vmeanscal;
166 for (
int i = 0; i < ttree->GetEntries(); ++i) {
169 if (dedx <= 0 || injtime < 0 || injring < 0)
continue;
171 double corrcetion =
getCorrection(injring, injtime, vmeanscal);
172 double old_cor =
m_DBInjectTime->getCorrection(
"mean", injring, injtime);
174 dedx = (dedx * old_cor) / corrcetion;
180 if (injring > 0.5) wr = 1;
183 unsigned int tb = htimes->GetXaxis()->FindBin(injtime);
190 double chi = (dedx - predmean) / predsigma;
191 hdedx_corr[wr][tb]->Fill(dedx);
192 hchi[wr][tb]->Fill(chi);
196 std::map<int, std::vector<double>> vmeans_chi, vresos_chi;
200 std::map<int, std::vector<double>> vresoscorr;
204 std::map<int, std::vector<double>> vresoscal;
205 std::array<double, numdedx::nrings> scale_reso;
209 std::map<int, std::vector<double>> vmeans_corr, vresos_corr;
217 for (
int ir = 0; ir < 2; ir++) {
259 B2INFO(
"dE/dx Injection time calibration done");
269 B2INFO(
"Saving calibration for: " <<
m_suffix <<
"");
275 std::map<
int, std::vector<double>>& vmeans, std::map<
int, std::vector<double>>& vresos)
277 for (
unsigned int ir = 0; ir <
c_rings; ir++) {
279 for (
unsigned int it = 0; it <
m_tbins; it++) {
280 double mean = 1.00, meanerr = 0.0;
281 double reso = 1.00, resoerr = 0.0;
282 if (hvar[ir][it]->Integral() > 250) {
285 if (status != fitOK) {
286 mean = hvar[ir][it]->GetMean();
287 hvar[ir][it]->SetTitle(Form(
"%s, (%d)", hvar[ir][it]->GetTitle(), status));
289 mean = hvar[ir][it]->GetFunction(
"gaus")->GetParameter(1);
290 meanerr = hvar[ir][it]->GetFunction(
"gaus")->GetParError(1);
291 reso = hvar[ir][it]->GetFunction(
"gaus")->GetParameter(2);
292 resoerr = hvar[ir][it]->GetFunction(
"gaus")->GetParError(2);
293 std::string title = Form(
"#mu_{fit}: %0.03f, #sigma_{fit}: %0.03f", mean, reso);
294 hvar[ir][it]->SetTitle(Form(
"%s, %s", hvar[ir][it]->GetTitle(), title.data()));
298 vmeans[ir * 2].push_back(mean);
299 vresos[ir * 2].push_back(reso);
300 vmeans[ir * 2 + 1].push_back(meanerr);
301 vresos[ir * 2 + 1].push_back(resoerr);
309 double histmean = temphist->GetMean();
310 double histrms = temphist->GetRMS();
311 temphist->GetXaxis()->SetRangeUser(histmean - 5.0 * histrms, histmean + 5.0 * histrms);
313 int fs = temphist->Fit(
"gaus",
"Q0");
315 B2INFO(Form(
"\tFit (round 1) for hist (%s) failed (status = %d)", temphist->GetName(), fs));
319 double mean = temphist->GetFunction(
"gaus")->GetParameter(1);
320 double width = temphist->GetFunction(
"gaus")->GetParameter(2);
321 temphist->GetXaxis()->SetRangeUser(mean - 5.0 * width, mean + 5.0 * width);
322 fs = temphist->Fit(
"gaus",
"QR",
"", mean -
m_sigmaR * width, mean +
m_sigmaR * width);
324 B2INFO(Form(
"\tFit (round 2) for hist (%s) failed (status = %d)", temphist->GetName(), fs));
328 temphist->GetXaxis()->SetRangeUser(mean - 5.0 * width, mean + 5.0 * width);
329 B2INFO(Form(
"\tFit for hist (%s) successful (status = %d)", temphist->GetName(), fs));
340 if (cruns == 0)B2INFO(
"CDCDedxInjectTimeAlgorithm: start exp " << expRun.first <<
" and run " << expRun.second <<
"");
345 int estart = erStart.first;
346 int rstart = erStart.second;
349 int rend = erEnd.second;
354 m_suffix = Form(
"e%dr%dto%d_nruns%d", estart, rstart, rend, cruns);
355 B2INFO(
"\t+ run = " << rend <<
", m_suffix = " <<
m_suffix <<
"");
357 m_suffix = Form(
"e%dr%d", estart, rstart);
358 B2INFO(
"tool run = " << estart <<
", exp = " << estart <<
", m_suffix = " <<
m_suffix <<
"");
370 for (
int ib = 0; ib < 69; ib++) {
371 fixedges[ib] = ib * 0.5 * 1e3;
372 if (ib > 40 && ib <= 60) fixedges[ib] = fixedges[ib - 1] + 1.0 * 1e3;
373 else if (ib > 60 && ib <= 64) fixedges[ib] = fixedges[ib - 1] + 10.0 * 1e3;
374 else if (ib > 64 && ib <= 65) fixedges[ib] = fixedges[ib - 1] + 420.0 * 1e3;
375 else if (ib > 65 && ib <= 66) fixedges[ib] = fixedges[ib - 1] + 500.0 * 1e3;
376 else if (ib > 66) fixedges[ib] = fixedges[ib - 1] + 2e6;
380 for (
unsigned int ib = 0; ib <
m_vtedges.size(); ib++)
388 for (
unsigned int ir = 0; ir <
c_rings; ir++) {
390 for (
unsigned int it = 0; it <
m_tbins; it++) {
392 std::string title = Form(
"%s(%s), time(%s)",
m_suffix.data(),
m_sring[ir].data(), label.data());
393 std::string hname = Form(
"h%s_%s_%s_t%d", var.data(),
m_sring[ir].data(),
m_suffix.data(), it);
396 else htemp[ir][it] =
new TH1D(hname.data(),
"", 50,
m_tedges[it],
m_tedges[it + 1]);
397 htemp[ir][it]->SetTitle(Form(
"%s;%s;entries", title.data(), var.data()));
406 const int nt[tzoom] = {50, 500, 1000};
407 double tzedges[tzoom + 1] = {0, 2.0e3, 0.1e6, 20e6};
408 std::string stname[tzoom] = {
"early",
"mid",
"later"};
409 std::string stlabel[tzoom] = {
"zoom <2ms",
"early time <= 100ms",
"later time >100ms"};
410 for (
unsigned int ir = 0; ir <
c_rings; ir++) {
411 for (
int wt = 0; wt < tzoom; wt++) {
412 std::string title = Form(
"inject time (%s), %s (%s)", stlabel[wt].data(),
m_sring[ir].data(),
m_suffix.data());
413 std::string hname = Form(
"htimezoom_%s_%s_%s",
m_sring[ir].data(), stname[wt].data(),
m_suffix.data());
414 htemp[ir][wt] =
new TH1D(Form(
"%s", hname.data()),
"", nt[wt], tzedges[wt], tzedges[wt + 1]);
415 htemp[ir][wt]->SetTitle(Form(
"%s;injection time(#mu-sec);entries", title.data()));
423 for (
unsigned int ir = 0; ir <
c_rings; ir++) {
424 for (
unsigned int it = 3; it <
m_tbins; it++) {
436 std::map<
int, std::vector<double>>& var,
437 std::map<
int, std::vector<double>>& time, TH1D*& htimes)
442 for (
int ir = 0; ir < 2; ir++) {
444 for (
int ix = varcorr[ir * 2].size(); ix -- > 0;) {
446 double var_thisbin = 1.0;
447 var_thisbin = var[ir * 2].at(ix);
449 double atime_thisbin = time[ir * 2].at(ix);
450 double ctime_thisbin = htimes->GetBinCenter(ix + 1);
452 if (atime_thisbin > 0 && atime_thisbin < 4e4 * 0.99) {
454 double var_nextbin = 1.0;
455 var_nextbin = var[ir * 2].at(ix + 1);
456 double var_diff = var_nextbin - var_thisbin;
458 double atime_nextbin = time[ir * 2].at(ix + 1);
459 double atime_diff = atime_nextbin - atime_thisbin;
461 double slope = (atime_diff > 0) ? var_diff / atime_diff : -1.0;
464 if (var_diff > 0 && slope > 0)varcorr[ir * 2].at(ix) = var_thisbin + (ctime_thisbin - atime_thisbin) * (slope);
465 printf(
"\t %s ix = %d, center = %0.2f(%0.3f), var = %0.5f(%0.5f) \n",
m_sring[ir].data(), ix, ctime_thisbin, atime_thisbin,
466 var_thisbin, varcorr[ir * 2].at(ix));
468 printf(
"\t %s --> ix = %d, center = %0.2f(%0.3f), var = %0.5f(%0.5f) \n",
m_sring[ir].data(), ix, ctime_thisbin, atime_thisbin,
469 var_thisbin, varcorr[ir * 2].at(ix));
477 std::map<
int, std::vector<double>>& var,
478 std::map<
int, std::vector<double>>& varscal, std::string svar)
482 B2INFO(
"CDCDedxInjectTimeAlgorithm: normalising constants with plateau");
483 for (
unsigned int ir = 0; ir <
c_rings; ir++) {
485 unsigned int msize = varscal[ir * 2].size();
488 for (
unsigned int im = 0; im < msize; im++) {
490 double mean = varscal[ir * 2].at(im);
491 if (time > 4e4 && mean > 0) {
496 if (countsum > 0 && scale[ir] > 0) {
497 scale[ir] /= countsum;
498 for (
unsigned int im = 0; im < msize; im++) {
499 varscal[ir * 2].at(im) /= scale[ir];
505 bool incomp_bin =
false;
506 std::vector<std::vector<double>> oldvectors;
508 int vsize = oldvectors.size();
509 if (vsize != 6) incomp_bin =
true;
511 for (
int iv = 0; iv < 2; iv++) {
512 if (oldvectors[iv * 3 + 1].size() != varscal[iv * 2].size()) incomp_bin =
true;
516 B2INFO(
"CDCDedxInjectTimeAlgorithm: started merging relative constants");
517 for (
int ir = 0; ir < 2; ir++) {
518 unsigned int msize = varscal[ir * 2].size();
519 for (
unsigned int im = 0; im < msize; im++) {
520 double relvalue = varscal[ir * 2].at(im);
521 double oldvalue = oldvectors[ir * 3 + 1].at(im);
522 double merged = oldvalue * relvalue;
523 printf(
"%s: rel %0.03f, old %0.03f, merged %0.03f\n",
m_suffix.data(), relvalue, oldvalue, merged);
524 varscal[ir * 2].at(im) *= oldvectors[ir * 3 + 1].at(im) ;
527 }
else B2ERROR(
"CDCDedxInjectTimeAlgorithm: found incompatible bins for merging");
528 }
else B2INFO(
"CDCDedxInjectTimeAlgorithm: saving final (abs) calibration");
534 TCanvas cfit(
"cfit",
"cfit", 1000, 500);
536 std::stringstream psname_fit;
537 psname_fit << Form(
"%s_%s_%s.pdf[",
m_prefix.data(), var.data(),
m_suffix.data());
538 cfit.Print(psname_fit.str().c_str());
540 psname_fit << Form(
"%s_%s_%s.pdf",
m_prefix.data(), var.data(),
m_suffix.data());
541 for (
unsigned int it = 0; it <
m_tbins; it++) {
542 for (
unsigned int ir = 0; ir <
c_rings; ir++) {
544 hvar[ir][it]->SetFillColorAlpha(ir + 5, 0.25);
545 hvar[ir][it]->Draw();
547 cfit.Print(psname_fit.str().c_str());
550 psname_fit << Form(
"%s_%s_%s.pdf]",
m_prefix.data(), var.data(),
m_suffix.data());
551 cfit.Print(psname_fit.str().c_str());
558 TCanvas cestat(
"cestat",
"cestat", 1000, 500);
559 cestat.SetBatch(kTRUE);
565 hestats->SetName(Form(
"hestats_%s",
m_suffix.data()));
566 hestats->SetStats(0);
567 hestats->Draw(
"hist text");
572 htstats->SetName(Form(
"htstats_%s",
m_suffix.data()));
573 htstats->SetStats(0);
574 htstats->Draw(
"hist text");
576 cestat.Print(Form(
"%s_eventstat_%s.pdf",
m_prefix.data(),
m_suffix.data()));
582 TCanvas ctzoom(
"ctzoom",
"ctzoom", 1500, 450);
583 ctzoom.SetBatch(kTRUE);
585 for (
int wt = 0; wt < 3; wt++) {
587 if (wt == 2) gPad->SetLogy();
588 for (
unsigned int ir = 0; ir <
c_rings; ir++) {
589 hvar[ir][wt]->SetStats(0);
590 hvar[ir][wt]->SetFillColorAlpha(5 + ir, 0.20);
592 double max1 = hvar[ir][wt]->GetMaximum();
593 double max2 = hvar[
c_rings - 1][wt]->GetMaximum();
594 if (max2 > max1) hvar[ir][wt]->SetMaximum(max2 * 1.05);
595 hvar[ir][wt]->Draw(
"");
596 }
else hvar[ir][wt]->Draw(
"same");
599 ctzoom.Print(Form(
"%s_timezoom_%s.pdf]",
m_prefix.data(),
m_suffix.data()));
604 std::map<
int, std::vector<double>>& vresos, std::map<
int, std::vector<double>>& corr, std::string svar)
606 std::string sname[3] = {
"mean",
"reso"};
607 const int lcolors[
c_rings] = {2, 4};
613 for (
int ic = 0; ic < 2; ic++) {
614 cconst[ic] =
new TCanvas(Form(
"c%sconst", sname[ic].data()),
"", 900, 500);
615 cconst[ic]->SetGridy(1);
618 TLegend* mleg =
new TLegend(0.50, 0.54, 0.80, 0.75, NULL,
"brNDC");
619 mleg->SetBorderSize(0);
620 mleg->SetFillStyle(0);
622 for (
unsigned int ir = 0; ir <
c_rings; ir++) {
624 std::string mtitle = Form(
"#mu(dedx), relative const. compare, (%s)",
m_suffix.data());
627 std::string rtitle = Form(
"#sigma(#chi), relative const. compare, (%s)",
m_suffix.data());
632 hmean[ir]->SetTitle(Form(
"%s;injection time(#mu-second);#mu (#chi-fit)", rtitle.data()));
633 hreso[ir]->SetTitle(Form(
"%s;injection time(#mu-second);#sigma (#chi-fit)", rtitle.data()));
634 hcorr[ir]->SetTitle(Form(
"%s;injection time(#mu-second);#sigma (#chi-fit bin-bais-corr)", rtitle.data()));
637 hmean[ir]->SetTitle(Form(
"%s;injection time(#mu-second);#mu (dedx-fit)", mtitle.data()));
638 hreso[ir]->SetTitle(Form(
"%s;injection time(#mu-second);#sigma (dedx-fit)", mtitle.data()));
639 hcorr[ir]->SetTitle(Form(
"%s;injection time(#mu-second);#mu (dedx-fit, bin-bais-corr)", mtitle.data()));
643 for (
unsigned int it = 0; it <
m_tbins; it++) {
647 hmean[ir]->SetBinContent(it + 1, vmeans[ir * 2].at(it));
648 hmean[ir]->SetBinError(it + 1, vmeans[ir * 2 + 1].at(it));
649 hmean[ir]->GetXaxis()->SetBinLabel(it + 1, label.data());
651 hreso[ir]->SetBinContent(it + 1, vresos[ir * 2].at(it));
652 hreso[ir]->SetBinError(it + 1, vresos[ir * 2 + 1].at(it));
653 hreso[ir]->GetXaxis()->SetBinLabel(it + 1, label.data());
655 hcorr[ir]->SetBinContent(it + 1, corr[ir * 2].at(it));
656 hcorr[ir]->SetBinError(it + 1, corr[ir * 2 + 1].at(it));
657 hcorr[ir]->GetXaxis()->SetBinLabel(it + 1, label.data());
660 mleg->AddEntry(hmean[ir], Form(
"%s",
m_sring[ir].data()),
"lep");
662 if (svar ==
"chi")
setHistStyle(hmean[ir], lcolors[ir], ir + 24, -0.60, 0.60);
663 else if (svar ==
"dedx")
setHistStyle(hmean[ir], lcolors[ir], ir + 24, 0.60, 1.10);
664 else setHistStyle(hmean[ir], lcolors[ir], ir + 24, 0.9, 1.10);
666 if (ir == 0)hmean[ir]->Draw(
"");
667 else hmean[ir]->Draw(
"same");
668 if (svar ==
"dedx") {
669 mleg->AddEntry(hcorr[ir], Form(
"%s (bin-bias-corr)",
m_sring[ir].data()),
"lep");
670 setHistStyle(hcorr[ir], lcolors[ir], ir + 20, 0.60, 1.10);
671 hcorr[ir]->Draw(
"same");
673 if (ir == 1)mleg->Draw(
"same");
676 if (svar ==
"chi")
setHistStyle(hreso[ir], lcolors[ir], ir + 24, 0.5, 1.50);
677 else setHistStyle(hreso[ir], lcolors[ir], ir + 24, 0.0, 0.15);
678 if (ir == 0)hreso[ir]->Draw(
"");
679 else hreso[ir]->Draw(
"same");
681 mleg->AddEntry(hcorr[ir], Form(
"%s (bin-bias-corr)",
m_sring[ir].data()),
"lep");
682 setHistStyle(hcorr[ir], lcolors[ir], ir + 20, 0.5, 1.50);
683 hcorr[ir]->Draw(
"same");
685 if (ir == 1)mleg->Draw(
"same");
688 for (
int ic = 0; ic < 2; ic++) {
689 cconst[ic]->SaveAs(Form(
"%s_relconst_%s_%s_%s.pdf",
m_prefix.data(), sname[ic].data(), svar.data(),
m_suffix.data()));
690 cconst[ic]->SaveAs(Form(
"%s_relconst_%s_%s_%s.root",
m_prefix.data(), sname[ic].data(), svar.data(),
m_suffix.data()));
693 for (
int ic = 0; ic < 2; ic++) {
704 const int lcolors[
c_rings] = {2, 4};
708 TCanvas* cconst =
new TCanvas(
"ctimestatconst",
"", 900, 500);
711 TLegend* rleg =
new TLegend(0.40, 0.60, 0.80, 0.72, NULL,
"brNDC");
712 rleg->SetBorderSize(0);
713 rleg->SetFillStyle(0);
715 for (
unsigned int ir = 0; ir <
c_rings; ir++) {
717 std::string title = Form(
"injection time, her-ler comparison, (%s)",
m_suffix.data());
719 htimestat[ir]->SetTitle(Form(
"%s;injection time(#mu-second);norm. entries", title.data()));
721 for (
unsigned int it = 0; it <
m_tbins; it++) {
723 htimestat[ir]->SetBinContent(it + 1, htime[ir][it]->Integral());
724 htimestat[ir]->SetBinError(it + 1, 0);
725 htimestat[ir]->GetXaxis()->SetBinLabel(it + 1, label.data());
729 double norm = htimestat[ir]->GetMaximum();
730 rleg->AddEntry(htimestat[ir], Form(
"%s (scaled with %0.02f)",
m_sring[ir].data(), norm),
"lep");
731 htimestat[ir]->Scale(1.0 / norm);
732 setHistStyle(htimestat[ir], lcolors[ir], ir + 20, 0.0, 1.10);
733 htimestat[ir]->SetFillColorAlpha(lcolors[ir], 0.30);
734 if (ir == 0) htimestat[ir]->Draw(
"hist");
735 else htimestat[ir]->Draw(
"hist same");
736 if (ir == 1)rleg->Draw(
"same");
739 cconst->SaveAs(Form(
"%s_relconst_timestat_%s.pdf",
m_prefix.data(),
m_suffix.data()));
740 cconst->SaveAs(Form(
"%s_relconst_timestat_%s.root",
m_prefix.data(),
m_suffix.data()));
746 std::map<
int, std::vector<double>>& vresos,
747 std::array<double, numdedx::nrings>& scale, std::array<double, numdedx::nrings>& scale_reso)
750 std::vector<std::vector<double>> oldvectors;
753 const int c_type = 2;
754 std::string sname[
c_rings] = {
"mean",
"reso"};
755 std::string stype[c_type] = {
"new",
"old"};
756 const int lcolors[
c_rings] = {2, 4};
757 const int lmarker[c_type] = {20, 24};
762 for (
int ic = 0; ic < 2; ic++) {
763 cconst[ic] =
new TCanvas(Form(
"c%sconst", sname[ic].data()),
"", 900, 500);
764 cconst[ic]->SetGridy(1);
767 TLegend* mleg =
new TLegend(0.50, 0.54, 0.80, 0.75, NULL,
"brNDC");
768 mleg->SetBorderSize(0);
769 mleg->SetFillStyle(0);
771 TLegend* rleg =
new TLegend(0.50, 0.54, 0.80, 0.75, NULL,
"brNDC");
772 rleg->SetBorderSize(0);
773 rleg->SetFillStyle(0);
775 for (
unsigned int ip = 0; ip < c_type; ip++) {
777 for (
unsigned int ir = 0; ir <
c_rings; ir++) {
779 std::string hname = Form(
"hfmean_%s_%s_%s",
m_sring[ir].data(), stype[ip].data(),
m_suffix.data());
781 std::string title = Form(
"#mu(dedx), final-mean-compare (%s)",
m_suffix.data());
782 hmean[ir][ip]->SetTitle(Form(
"%s;injection time(#mu-second);#mu (dedx-fit)", title.data()));
784 hname = Form(
"hfreso_%s_%s_%s",
m_sring[ir].data(), stype[ip].data(),
m_suffix.data());
786 title = Form(
"#sigma(#chi), final-reso-compare (%s)",
m_suffix.data());
787 hreso[ir][ip]->SetTitle(Form(
"%s;injection time(#mu-second);#sigma (#chi-fit)", title.data()));
789 for (
unsigned int it = 0; it <
m_tbins; it++) {
792 double mean = 0.0, reso = 0.0;
799 mean = oldvectors[ir * 3 + 1].at(it);
800 reso = oldvectors[ir * 3 + 2].at(it);
805 hmean[ir][ip]->SetBinContent(it + 1, mean);
806 hmean[ir][ip]->SetBinError(it + 1, vmeans[ir * 2 + 1].at(it));
807 hmean[ir][ip]->GetXaxis()->SetBinLabel(it + 1, label.data());
809 hreso[ir][ip]->SetBinContent(it + 1, reso);
810 hreso[ir][ip]->SetBinError(it + 1, vresos[ir * 2 + 1].at(it));
811 hreso[ir][ip]->GetXaxis()->SetBinLabel(it + 1, label.data());
815 if (ip == 1)mleg->AddEntry(hmean[ir][ip], Form(
"%s, %s",
m_sring[ir].data(), stype[ip].data()),
"lep");
816 else mleg->AddEntry(hmean[ir][ip], Form(
"%s, %s (scaled by %0.03f)",
m_sring[ir].data(), stype[ip].data(), scale[ir]),
"lep");
817 setHistStyle(hmean[ir][ip], lcolors[ir], lmarker[ip] + ir * 2, 0.60, 1.05);
818 if (ir == 0 && ip == 0)hmean[ir][ip]->Draw(
"");
819 else hmean[ir][ip]->Draw(
"same");
820 if (ir == 1 && ip == 1)mleg->Draw(
"same");
823 if (ip == 1)rleg->AddEntry(hreso[ir][ip], Form(
"%s, %s",
m_sring[ir].data(), stype[ip].data()),
"lep");
824 else rleg->AddEntry(hreso[ir][ip], Form(
"%s, %s (scaled by %0.03f)",
m_sring[ir].data(), stype[ip].data(), scale_reso[ir]),
"lep");
825 setHistStyle(hreso[ir][ip], lcolors[ir], lmarker[ip] + ir * 2, 0.6, 1.9);
826 if (ir == 0 && ip == 0)hreso[ir][ip]->Draw(
"");
827 else hreso[ir][ip]->Draw(
"same");
828 if (ir == 1 && ip == 1)rleg->Draw(
"same");
832 for (
int ic = 0; ic < 2; ic++) {
833 cconst[ic]->SaveAs(Form(
"%s_finalconst_%s_%s.pdf",
m_prefix.data(), sname[ic].data(),
m_suffix.data()));
834 cconst[ic]->SaveAs(Form(
"%s_finalconst_%s_%s.root",
m_prefix.data(), sname[ic].data(),
m_suffix.data()));
838 for (
unsigned int ir = 0; ir <
c_rings; ir++) {
839 for (
unsigned int ip = 0; ip < c_type; ip++) {
840 delete hmean[ir][ip];
841 delete hreso[ir][ip];
848 std::map<
int, std::vector<double>>& vmeans)
851 unsigned int iv = ring * 2;
855 std::vector<unsigned int> tedges(sizet);
858 if (time >= 5e6) time = 5e6 - 10;
866 int thisbin = it, nextbin = it;
867 if (center != time && it > 0) {
872 if (it < sizet - 2)nextbin = it + 1;
873 else thisbin = it - 1;
877 double diff = vmeans[iv].at(2) - vmeans[iv].at(1) ;
880 if (it == 1) nextbin = it;
881 else nextbin = it + 1;
891 double thisdedx = vmeans[iv].at(thisbin);
892 double nextdedx = vmeans[iv].at(nextbin);
897 double newdedx = vmeans[iv].at(it);
898 if (thisbin != nextbin)
899 newdedx = thisdedx + ((nextdedx - thisdedx) / (nexttime - thistime)) * (time - thistime);
bool m_isMerge
merge payload when rel constant
void plotBinLevelDist(std::array< std::vector< TH1D * >, numdedx::nrings > &hvar, std::string var)
function to draw dedx, chi and time dist.
std::vector< double > m_vtlocaledges
internal time vector
double m_sigmaR
fit dedx dist in sigma range
void correctBinBias(std::map< int, std::vector< double > > &varcorr, std::map< int, std::vector< double > > &var, std::map< int, std::vector< double > > &time, TH1D *&htimes)
function to correct dedx mean/reso and return corrected vector map
void defineHisto(std::array< std::vector< TH1D * >, numdedx::nrings > &htemp, std::string var)
function to define histograms for dedx and time dist.
double getCorrection(unsigned int ring, unsigned int time, std::map< int, std::vector< double > > &vmeans)
function to get the correction factor of mean
double m_chiMax
max range of chi
std::string m_prefix
string prefix for plot names
void getExpRunInfo()
function to get exp/run information (payload object, plotting)
void setHistStyle(TH1D *&htemp, const int ic, const int is, const double min, const double max)
function to set histogram cosmetics
void plotFinalConstants(std::map< int, std::vector< double > > &vmeans, std::map< int, std::vector< double > > &vresos, std::array< double, numdedx::nrings > &scale, std::array< double, numdedx::nrings > &scale_reso)
function to final constant from merging or abs fits
int m_countR
a hack for running functions once
double * m_tedges
internal time array (copy of vtlocaledges)
double m_chiMin
min range of chi
bool m_ismakePlots
produce plots for monitoring
std::string m_suffix
string suffix for object names
int m_chiBins
bins for chi histogram
static const int c_rings
injection ring constants
int m_dedxBins
bins for dedx histogram
void createPayload(std::array< double, numdedx::nrings > &scale, std::map< int, std::vector< double > > &vmeans, std::map< int, std::vector< double > > &varscal, std::string svar)
function to store payloads after full calibration
void plotInjectionTime(std::array< std::array< TH1D *, 3 >, numdedx::nrings > &hvar)
function to injection time distributions (HER/LER in three bins)
void plotTimeStat(std::array< std::vector< TH1D * >, numdedx::nrings > &htime)
function to draw time stats
void fitGaussianWRange(TH1D *&temphist, fstatus &status)
function to perform gauss fit for input histogram
void defineTimeBins()
function to set/reset time bins
void plotRelConstants(std::map< int, std::vector< double > > &vmeans, std::map< int, std::vector< double > > &vresos, std::map< int, std::vector< double > > &corr, std::string svar)
function to relative constant from dedx fit mean and chi fit reso
void defineTimeHisto(std::array< std::array< TH1D *, 3 >, numdedx::nrings > &htemp)
function to define injection time bins histograms (monitoring only)
void deleteHisto(std::array< std::vector< TH1D * >, numdedx::nrings > &htemp)
function to delete histograms for dedx and time dist.
void deleteTimeHisto(std::array< std::array< TH1D *, 3 >, numdedx::nrings > &htemp)
function to define injection time bins histograms (monitoring only)
void plotEventStats()
function to draw event/track statistics plots
virtual EResult calibrate() override
CDC dE/dx Injection time algorithm.
void checkStatistics(std::array< std::vector< TH1D * >, numdedx::nrings > &hvar)
check statistics for obtaining calibration const.
std::vector< double > m_vtedges
external time vector
std::vector< std::vector< double > > m_vinjPayload
vector to store payload values
double m_dedxMax
max range of dedx
bool m_isminStat
flag to merge runs for statistics thershold
CDCDedxInjectTimeAlgorithm()
Constructor: Sets the description, the properties and the parameters of the algorithm.
std::array< std::string, numdedx::nrings > m_sring
injection ring name
DBObjPtr< CDCDedxInjectionTime > m_DBInjectTime
Injection time DB object.
double m_dedxMin
min range of dedx
unsigned int m_tbins
internal time bins
int m_thersE
min tracks to start calibration
void getMeanReso(std::array< std::vector< TH1D * >, numdedx::nrings > &hvar, std::map< int, std::vector< double > > &vmeans, std::map< int, std::vector< double > > &vresos)
function to get mean and reso of histogram
std::string getTimeBinLabel(const double &tedges, const int &it)
function to return time label for histograms labeling
dE/dx injection time calibration constants
Class to hold the prediction of mean as a function of beta-gamma (bg)
double getMean(double bg)
Return the predicted mean value as a function of beta-gamma (bg)
void setParameters(std::string infile)
set the parameters from file
Class to hold the prediction of resolution depending dE/dx, nhit, and cos(theta)
double ionzPrediction(double dedx)
Return sigma from the ionization parameterization.
double cosPrediction(double cos)
Return sigma from the cos parameterization.
double nhitPrediction(double nhit)
Return sigma from the nhit parameterization.
void setParameters(std::string infile)
set the parameters from file
void saveCalibration(TClonesArray *data, const std::string &name)
Store DBArray payload with given name with default IOV.
void updateDBObjPtrs(const unsigned int event, const int run, const int experiment)
Updates any DBObjPtrs by calling update(event) for DBStore.
void setDescription(const std::string &description)
Set algorithm description (in constructor)
const std::vector< Calibration::ExpRun > & getRunList() const
Get the list of runs for which calibration is called.
EResult
The result of calibration.
@ c_OK
Finished successfully =0 in Python.
@ c_NotEnoughData
Needs more data =2 in Python.
CalibrationAlgorithm(const std::string &collectorModuleName)
Constructor - sets the prefix for collected objects (won't be accesses until execute(....
static const ChargedStable electron
electron particle
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...
Abstract base class for different kinds of events.