52 double dedx, costh;
int charge;
53 ttree->SetBranchAddress(
"dedx", &dedx);
54 ttree->SetBranchAddress(
"costh", &costh);
55 ttree->SetBranchAddress(
"charge", &charge);
58 vector<TH1D*> hdedx_negi, hdedx_posi;
63 for (
unsigned int i = 0; i <
m_npBins; ++i) {
65 double mincos = i * bwposi +
m_posMin;
66 double maxcos = mincos + bwposi;
67 string title = Form(
"costh: %0.04f, %0.04f(%s)", mincos, maxcos,
m_suffix.data());
69 hdedx_posi[i]->SetTitle(Form(
"%s;dedx;entries", title.data()));
72 maxcos = mincos + bwnegi;
73 title = Form(
"costh: %0.04f, %0.04f(%s)", mincos, maxcos,
m_suffix.data());
75 hdedx_negi[i]->SetTitle(Form(
"%s;dedx;entries", title.data()));
79 for (
int i = 0; i < ttree->GetEntries(); ++i) {
84 if (dedx <= 0 || charge == 0)
continue;
85 if (costh > -0.850 && costh < 0.950)
continue;
88 icosbin = int((costh -
m_posMin) / bwposi);
89 hdedx_posi[icosbin]->Fill(dedx);
90 }
else if (costh < 0) {
91 icosbin = int((costh -
m_negMin) / bwnegi);
92 hdedx_negi[icosbin]->Fill(dedx);
96 map<int, vector<double>> vneg_fitpars;
97 map<int, vector<double>> vpos_fitpars;
99 vector<double> vneg_const, vpos_const;
100 vector<vector<double>> vfinal_const;
102 for (
unsigned int i = 0; i <
m_npBins; ++i) {
108 if (status != FitOK) {
109 vneg_fitpars[0].push_back(1.0);
110 vneg_fitpars[1].push_back(0.0);
111 vneg_fitpars[2].push_back(0.0);
112 vneg_fitpars[3].push_back(0.0);
113 hdedx_negi[i]->SetTitle(Form(
"%s, Fit(%d)", hdedx_negi[i]->GetTitle(), status));
115 vneg_fitpars[0].push_back(hdedx_negi[i]->GetFunction(
"gaus")->GetParameter(1));
116 vneg_fitpars[1].push_back(hdedx_negi[i]->GetFunction(
"gaus")->GetParError(1));
117 vneg_fitpars[2].push_back(hdedx_negi[i]->GetFunction(
"gaus")->GetParameter(2));
118 vneg_fitpars[3].push_back(hdedx_negi[i]->GetFunction(
"gaus")->GetParError(2));
121 vneg_const.push_back(vneg_fitpars[0][i]);
125 if (status != FitOK) {
126 vpos_fitpars[0].push_back(1.0);
127 vpos_fitpars[1].push_back(0.0);
128 vpos_fitpars[2].push_back(0.0);
129 vpos_fitpars[3].push_back(0.0);
130 hdedx_posi[i]->SetTitle(Form(
"%s, Fit(%d)", hdedx_posi[i]->GetTitle(), status));
132 vpos_fitpars[0].push_back(hdedx_posi[i]->GetFunction(
"gaus")->GetParameter(1));
133 vpos_fitpars[1].push_back(hdedx_posi[i]->GetFunction(
"gaus")->GetParError(1));
134 vpos_fitpars[2].push_back(hdedx_posi[i]->GetFunction(
"gaus")->GetParameter(2));
135 vpos_fitpars[3].push_back(hdedx_posi[i]->GetFunction(
"gaus")->GetParError(2));
138 vpos_const.push_back(vpos_fitpars[0][i]);
141 vfinal_const.push_back(vneg_const);
142 vfinal_const.push_back(vpos_const);
149 plotHist(hdedx_posi, vpos_fitpars,
"pos");
150 plotHist(hdedx_negi, vneg_fitpars,
"neg");
194 B2FATAL(
"CDCDedxCosEdgeAlgorithm: Can't merge paylaods with different size");
196 for (
unsigned int ibin = 0; ibin <
m_npBins; ibin++) {
200 double relg = vfinalconst[0].at(ibin);
201 double newg = prevg * relg;
202 B2INFO(
"CosEdge Const (<0), bin# " << ibin <<
", rel " << relg <<
", previous " << prevg <<
", merged " << newg);
203 vfinalconst[0].at(ibin) *= (double)
m_DBCosineCor->getMean(-1, ibin);
204 vfinalconst[0].at(ibin) /= (0.5 * (vfinalconst[0].at(
m_npBins - 1) + vfinalconst[0].at(
m_npBins - 2)));
208 relg = vfinalconst[1].at(ibin);
210 B2INFO(
"CosEdge Const (>0), bin# " << ibin <<
", rel " << relg <<
", previous " << prevg <<
", merged " << newg);
211 vfinalconst[1].at(ibin) *= (double)
m_DBCosineCor->getMean(1, ibin);
212 vfinalconst[1].at(ibin) /= (0.5 * (vfinalconst[1].at(0) + vfinalconst[1].at(1)));
216 B2INFO(
"CDCDedxCosineEdge calibration done");
224 if (temphist->Integral() < 500) {
225 B2INFO(Form(
"\t insufficient fit stats (%0.00f) for (%s)", temphist->Integral(), temphist->GetName()));
229 temphist->GetXaxis()->SetRange(temphist->FindFirstBinAbove(0, 1), temphist->FindLastBinAbove(0, 1));
230 int fs = temphist->Fit(
"gaus",
"QR");
232 B2INFO(Form(
"\tFit (round 1) for hist (%s) failed (status = %d)", temphist->GetName(), fs));
236 double fdEdxMean = temphist->GetFunction(
"gaus")->GetParameter(1);
237 double width = temphist->GetFunction(
"gaus")->GetParameter(2);
238 temphist->GetXaxis()->SetRangeUser(fdEdxMean - 5.0 * width, fdEdxMean + 5.0 * width);
239 fs = temphist->Fit(
"gaus",
"QR",
"", fdEdxMean -
m_sigLim * width, fdEdxMean +
m_sigLim * width);
241 B2INFO(Form(
"\tFit (round 2) for hist (%s) failed (status = %d)", temphist->GetName(), fs));
245 temphist->GetXaxis()->SetRangeUser(fdEdxMean - 5.0 * width, fdEdxMean + 5.0 * width);
246 B2INFO(Form(
"\tFit for hist (%s) successful (status = %d)", temphist->GetName(), fs));
257 TCanvas ctmp(
"ctmp",
"ctmp", 1200, 1200);
261 psname << Form(
"cdcdedx_cosedgecal_fits_%s_%s.pdf[", type.data(),
m_suffix.data());
262 ctmp.Print(psname.str().c_str());
264 psname << Form(
"cdcdedx_cosedgecal_fits_%s_%s.pdf", type.data(),
m_suffix.data());
266 for (
unsigned int i = 0; i <
m_npBins; ++i) {
269 hdedx[i]->SetStats(0);
270 hdedx[i]->SetFillColor(kAzure - 9);
271 hdedx[i]->DrawCopy();
273 TPaveText* pt =
new TPaveText(0.5, 0.73, 0.8, 0.89,
"NBNDC");
275 pt->AddText(Form(
"#mu_{fit}: %0.03f#pm%0.03f", vpars[0][i], vpars[1][i]));
276 pt->AddText(Form(
"#sigma_{fit}: %0.03f#pm%0.03f", vpars[2][i], vpars[3][i]));
279 if ((i + 1) % 20 == 0 || (i + 1) ==
m_npBins) {
280 ctmp.Print(psname.str().c_str());
288 psname << Form(
"cdcdedx_cosedgecal_fits_%s_%s.pdf]", type.data(),
m_suffix.data());
289 ctmp.Print(psname.str().c_str());
294 map<
int, vector<double>>& vpos_fitpars)
297 TCanvas cQa(
"cQa",
"cQa", 1200, 1200);
300 double min[2] = {0.85, 0.04};
301 double max[2] = {1.05, 0.3};
303 string vars[2] = {
"#mu_{fit}",
"#sigma_{fit}"};
304 string side[2] = {
"pcos",
"ncos"};
306 for (
int is = 0; is < 2; is++) {
314 for (
int iv = 0; iv < 2; iv++) {
316 string hname = Form(
"hpar_%s_%s_%s", vars[iv].data(), side[is].data(),
m_suffix.data());
318 TH1D htemp(Form(
"%s", hname.data()),
"",
m_npBins, minp, maxp);
319 htemp.SetTitle(Form(
"Constant (%s), cos#theta:(%0.02f, %0.02f);cos(#theta);const", vars[iv].data(), minp, maxp));
321 for (
unsigned int ib = 0; ib <
m_npBins; ib++) {
323 htemp.SetBinContent(ib + 1, vpos_fitpars[2 * iv][ib]);
324 htemp.SetBinError(ib + 1, vpos_fitpars[2 * iv + 1][ib]);
326 htemp.SetBinContent(ib + 1, vneg_fitpars[2 * iv][ib]);
327 htemp.SetBinError(ib + 1, vneg_fitpars[2 * iv + 1][ib]);
333 cQa.cd(2 * iv + 1 + is);
339 cQa.SaveAs(Form(
"cdcdedx_cosedgecal_relconst_%s.pdf",
m_suffix.data()));
340 cQa.SaveAs(Form(
"cdcdedx_cosedgecal_relconst_%s.root",
m_suffix.data()));
347 TCanvas ctmp_const(
"ctmp_const",
"ctmp_const", 900, 450);
348 ctmp_const.Divide(2, 1);
350 for (
int i = 0; i < 2; i++) {
358 TH1D holdconst(Form(
"holdconst%d_%s", i,
m_suffix.data()),
"",
m_npBins, min, max);
359 holdconst.SetTitle(Form(
"constant comparison, cos#theta:(%0.02f, %0.02f);cos(#theta);const", min, max));
361 TH1D hnewconst(Form(
"hnewconst%d_%s", i,
m_suffix.data()),
"",
m_npBins, min, max);
363 int iside = 2 * i - 1;
364 for (
int ibin = 0; ibin <
m_DBCosineCor->getSize(iside); ibin++) {
365 holdconst.SetBinContent(ibin + 1, (
double)
m_DBCosineCor->getMean(iside, ibin));
366 hnewconst.SetBinContent(ibin + 1, vfinalconst[i].at(ibin));
369 ctmp_const.cd(i + 1);
374 holdconst.DrawCopy(
"");
377 hnewconst.DrawCopy(
"same");
379 TPaveText* pt =
new TPaveText(0.47, 0.73, 0.77, 0.89,
"NBNDC");
382 TText* told = pt->AddText(
"old const");
383 told->SetTextColor(kBlack);
384 TText* tnew = pt->AddText(
"new const");
385 tnew->SetTextColor(kRed);
390 ctmp_const.SaveAs(Form(
"cdcdedx_cosedgecal_constcomp_%s.pdf",
m_suffix.data()));
391 ctmp_const.SaveAs(Form(
"cdcdedx_cosedgecal_constcomp_%s.root",
m_suffix.data()));